
#include "ObjectTemplate.h"

#define  PRINT_ANALYSIS_IMAGE  1

void ideal_observer::DrawGridImage(CImg<int>& imgOutput, const CImg<int>& imgContext, int IW, int IH, int CW, int CH, int H_Extent)
{
	if (imgContext.width != IW * CW) return;
	if (imgContext.height != IH * CH) return;
	int iiw = (IW + 1) * CW;
	int iih = (IH + 1) * CH;
	imgOutput.assign(iiw, iih + H_Extent, 1, 3);
	imgOutput.fill(255);

	int ix, iy;
	for ( iy = CH; iy < iih; ++iy) for ( ix = CW; ix < iiw; ++ix )
	{
		imgOutput.at(ix, iy, 0, 0) = imgContext.at(ix - CW, iy - CH, 0, (imgContext.dim>0 ? 0 : imgContext.dim-1));
		imgOutput.at(ix, iy, 0, 1) = imgContext.at(ix - CW, iy - CH, 0, (imgContext.dim>1 ? 1 : imgContext.dim-1));
		imgOutput.at(ix, iy, 0, 2) = imgContext.at(ix - CW, iy - CH, 0, (imgContext.dim>2 ? 2 : imgContext.dim-1));
	}
	//for ( iy = 0; iy <= IH; ++iy ) imgOutput.draw_line(0, iy * CH + CH - 1, iiw, iy * CH + CH - 1, CExternalTools::GetIntArrayRGB(192, 192, 192));
	//for ( ix = 0; ix <= IW; ++ix ) imgOutput.draw_line(ix * CW + CW - 1, 0, ix * CW + CW - 1, iih, CExternalTools::GetIntArrayRGB(192, 192, 192));

	int cc, cd, fs, ox, oy;
	double tfs;
	char strTemp[10001];

	cc = 1; cd = 10; while (IW>cd) { ++cc; cd *= 10; }
	tfs = 2.0 * CW / cc; if (tfs > CH) tfs = CH;
	fs = (int)(tfs * 0.6);
	ox = (CW-fs/2*cc) / 2;
	oy = (CH-fs) / 2;

	for ( iy = 1; iy <= IH; ++iy )
	{
		sprintf_s(strTemp, "%0*d", cc, iy - 1);
		CExternalTools::OutputTextToImage(imgOutput, strTemp, ox, oy + iy * CH, 255, 255, 255, 128, 128, 255, fs);
	}
	for ( ix = 1; ix <= IW; ++ix )
	{
		sprintf_s(strTemp, "%0*d", cc, ix - 1);
		CExternalTools::OutputTextToImage(imgOutput, strTemp, ox + ix * CW, oy, 255, 255, 255, 128, 128, 255, fs);
	}
	for ( iy = 0; iy < H_Extent; ++iy ) for ( ix = 0; ix < iiw; ++ix )
		imgOutput.at(ix, iih + iy, 0, 0) = imgOutput.at(ix, iih + iy, 0, 1) = imgOutput.at(ix, iih + iy, 0, 2) = 0;
}

// [*] CObjectEdglet

// [*] CTransformation

ideal_observer::CTransformation::CTransformation()
{
	Clear();
}

ideal_observer::CTransformation::CTransformation(int base_x, int base_y)
{
	myInit(base_x, base_y);
}

void ideal_observer::CTransformation::myInit(int base_x, int base_y)
{
	this->base_x = base_x;
	this->base_y = base_y;
}

void ideal_observer::CTransformation::Clear()
{
	base_x = base_y = 0;
}

bool ideal_observer::CTransformation::InputFromStream( istream& inStream )
{
	if (!(inStream>>base_x>>base_y)) return false;
	return true;
}

void ideal_observer::CTransformation::OutputToStream( ofstream& outStream ) const
{
	outStream<<"\t"<<base_x;
	outStream<<"\t"<<base_y;
	outStream<<endl;
}

// CTransformationSet

ideal_observer::CTransformationSet::CTransformationSet()
{
	Clear();
}

void ideal_observer::CTransformationSet::Clear()
{
	Size = 0;
	List.Clear();
	Weights.Clear();
}

void ideal_observer::CTransformationSet::PushBack(const CTransformation& NewTransformation, double NewWeight)
{
	++Size;
	List.PushBack(NewTransformation);
	Weights.PushBack(NewWeight);
}

void ideal_observer::CTransformationSet::NormalizeWeights()
{
	int i;
	double ts = 0.0;
	for ( i = 0; i < Size; ++i )
		ts += Weights.pElement[i];
	if (ts < EPS) return;
	for ( i = 0; i < Size; ++i )
		Weights.pElement[i] /= ts;
}

bool ideal_observer::CTransformationSet::InputFromFile( string strFN_SrcSet )
{
	int i;
	bool ret = true;

	ifstream inFile(strFN_SrcSet.c_str());
	ret = (inFile>>Size) && (Size>=0);
	Size = (Size<0 ? 0 : Size);
	List.Resize(Size);
	Weights.Resize(Size);
	for ( i = 0; ret && i<Size; ++i )
		ret = ret
			&& (inFile>>Weights.pElement[i])
			&& List.pElement[i].InputFromStream(inFile);
	inFile.close();
	inFile.clear();

	if (!ret) List.Clear();
	NormalizeWeights();
	return ret;
}

void ideal_observer::CTransformationSet::OutputToFile( string strFN_DstSet )
{
	int i;

	ofstream outFile(strFN_DstSet.c_str());
	outFile<<Size<<endl;
	for ( i = 0; i < Size; ++i )
	{
		outFile<<Weights.pElement[i]<<"\t";
		List.pElement[i].OutputToStream(outFile);
	}
	outFile.close();
	outFile.clear();
}

// [*] CFixation

ideal_observer::CFixation::CFixation()
{
	Clear();
}

ideal_observer::CFixation::CFixation(int ob_x, int ob_y)
{
	myInit(ob_x, ob_y);
}

void ideal_observer::CFixation::myInit(int ob_x, int ob_y)
{
	this->ob_x = ob_x;
	this->ob_y = ob_y;
}

void ideal_observer::CFixation::Clear()
{
	ob_x = ob_y = 0;
}

bool ideal_observer::CFixation::InputFromStream( istream& inStream )
{
	if (!(inStream>>ob_x>>ob_y)) return false;
	return true;
}

void ideal_observer::CFixation::OutputToStream( ofstream& outStream ) const
{
	outStream<<"\t"<<ob_x;
	outStream<<"\t"<<ob_y;
	outStream<<endl;
}

// CFixationSet

ideal_observer::CFixationSet::CFixationSet()
{
	Clear();
}

void ideal_observer::CFixationSet::Clear()
{
	Size = 0;
	List.Clear();
	//Weights.Clear();
}

void ideal_observer::CFixationSet::PushBack(const CFixation& NewFixation) //, double NewWeight)
{
	++Size;
	List.PushBack(NewFixation);
	//Weights.PushBack(NewWeight);
}

//void ideal_observer::CFixationSet::NormalizeWeights()
//{
//	int i;
//	double ts = 0.0;
//	for ( i = 0; i < Size; ++i )
//		ts += Weights.pElement[i];
//	if (ts < EPS) return;
//	for ( i = 0; i < Size; ++i )
//		Weights.pElement[i] /= ts;
//}

bool ideal_observer::CFixationSet::InputFromFile( string strFN_SrcSet )
{
	int i;
	bool ret = true;

	ifstream inFile(strFN_SrcSet.c_str());
	ret = (inFile>>Size) && (Size>=0);
	Size = (Size<0 ? 0 : Size);
	List.Resize(Size);
	//Weights.Resize(Size);
	for ( i = 0; ret && i<Size; ++i )
		ret = ret
			//&& (inFile>>Weights.pElement[i])
			&& List.pElement[i].InputFromStream(inFile);
	inFile.close();
	inFile.clear();

	if (!ret) List.Clear();
	//NormalizeWeights();
	return ret;
}

void ideal_observer::CFixationSet::OutputToFile( string strFN_DstSet )
{
	int i;

	ofstream outFile(strFN_DstSet.c_str());
	outFile<<Size<<endl;
	for ( i = 0; i < Size; ++i )
	{
		//outFile<<Weights.pElement[i]<<"\t";
		List.pElement[i].OutputToStream(outFile);
	}
	outFile.close();
	outFile.clear();
}

// [*] CObjectTemplate

void ideal_observer::CObjectTemplate::Clear()
{
	Amount_Edglet = 0;
	List_Edglet.Clear();
}

void ideal_observer::CObjectTemplate::myInit(const CSimpleTypeArray<int>& matSrcObserveResult, int IW, int IH)
{
	int ix, iy;
	CObjectEdglet toe;

	this->IW = IW;
	this->IH = IH;
	this->IWH = IW * IH;

	xLow = IW; yLow = IH;
	xHigh = -1; yHigh = -1;

	xCenter = yCenter = 0;
	this->Clear();
	for ( iy = 0; iy < IH; ++iy ) for ( ix = 0; ix < IW; ++ix )
	{
		toe.x = ix; toe.y = iy;
		toe.edglet.myInit(matSrcObserveResult.pElement[ix + iy * IW]);
		if (!(toe.edglet.fDoesExist())) continue;

		++Amount_Edglet;
		List_Edglet.PushBack(toe);
		xCenter += ix; yCenter += iy;
		if (ix < xLow) xLow = ix;
		if (ix > xHigh) xHigh = ix;
		if (iy < yLow) yLow = iy;
		if (iy > yHigh) yHigh = iy;
	}
	if (Amount_Edglet > 0)
	{
		xCenter = (xCenter + Amount_Edglet / 2) / Amount_Edglet;
		yCenter = (yCenter + Amount_Edglet / 2) / Amount_Edglet;
	}
}

void ideal_observer::CObjectTemplate::DoTransformation(
	CSimpleTypeArray<CObjectEdglet>& DstListOE,
	const CTransformation& _Transformation,
	bool fDeleteEmptyEdglets) const
{
	DstListOE.Clear();

	int i_oe, ix, iy;
	CObjectEdglet toe;

	if (!fDeleteEmptyEdglets)
	{
		DstListOE.Resize(IWH);
		for ( iy = 0; iy < IH; ++iy ) for ( ix = 0; ix < IW; ++ix )
		{
			CObjectEdglet& cur_oe = DstListOE.pElement[ix + iy * IW];
			cur_oe.x = ix;
			cur_oe.y = iy;
			cur_oe.edglet.Clear();
		}
	}

	for ( i_oe = 0; i_oe < Amount_Edglet; ++i_oe )
	{
		toe = List_Edglet.pElement[i_oe];
		toe.x += _Transformation.base_x;
		toe.y += _Transformation.base_y;
		if (toe.x<0 || toe.x>=IW || toe.y<0 || toe.y>=IH) continue;
		if (fDeleteEmptyEdglets)
			DstListOE.PushBack(toe);
		else
			DstListOE.pElement[toe.x + toe.y * IW] = toe;
	}
}

void ideal_observer::CObjectTemplate::DoTransformation( CSimpleTypeArray<int>& matD, const CTransformation& _Transformation ) const
{
	matD.Resize(IWH);
	matD.Fill(-1);
	int i_oe, ix, iy;
	for ( i_oe = 0; i_oe < Amount_Edglet; ++i_oe )
	{
		const CObjectEdglet& toe = List_Edglet.pElement[i_oe];
		ix = toe.x + _Transformation.base_x;
		iy = toe.y + _Transformation.base_y;
		if (ix<0 || ix>=IW || iy<0 || iy>=IH) continue;
		matD.pElement[ix + iy * IW] = toe.edglet.label;
	}
}

// [*] CDistribution

ideal_observer::CDistribution::CDistribution()
{
	myInit(0, 0);
}

ideal_observer::CDistribution::CDistribution( int IW, int IH )
{
	myInit(IW, IH);
}

void ideal_observer::CDistribution::myInit( int IW, int IH )
{
	IE = CEdglet::Amount_Edglet;
	this->IW = IW;
	this->IH = IH;
	this->IWH = IW * IH;
	this->Size = ((size_t)IE) * (size_t)IWH;
	ce_logP.Resize(Size);
	ce_logP.Fill(-CMyMath::my_log(IE));
	ce_realP.Resize(Size);
	ce_realP.Fill(1.0/IE);
}

void ideal_observer::CDistribution::ClearProb()
{
	ce_logP.Fill(-CMyMath::my_log(IE));
	ce_realP.Fill(1.0/IE);
}

void ideal_observer::CDistribution::DeepClear()
{
	myInit(0, 0);
	ce_logP.DeepClear();
	ce_realP.DeepClear();
}

void ideal_observer::CDistribution::OutputToImage( CImg<int>& imgOutput, int CW, int CH, int H_Extent, bool fDrawGrid ) const
{
	if (IWH < 1) return;

	int i;
	CRectBox<int> iBox;
	int iiw, iih;
	CImg<int> imgContext;
	iiw = IW * CW; iih = IH * CH;
	imgContext.assign(iiw, iih, 1, 1);
	imgContext.fill(255);
	const double* qd;
	iBox.w = CW; iBox.h = CH;
	for ( i = 0; i < IWH; ++i )
	{
		iBox.x = i % IW * CW;
		iBox.y = i / IW * CH;
		qd = this->pBase_realP(i);
		CEdglet::DrawCellDistributionToImage(imgContext.data, iiw, iih, iBox, qd);
	}

	if (fDrawGrid)
		DrawGridImage(imgOutput, imgContext, IW, IH, CW, CH, H_Extent);
	else
	{
		imgOutput.assign(iiw, iih + H_Extent, 1, 3);
		imgOutput.fill(0);
		imgOutput.draw_image(imgContext, 0, 0, 0, 0);
		imgOutput.draw_image(imgContext, 0, 0, 0, 1);
		imgOutput.draw_image(imgContext, 0, 0, 0, 2);
	}
}

double* ideal_observer::CDistribution::pBase_logP( int a_offset )
{
	if (a_offset<0 || a_offset>=IWH) return NULL;
	return ce_logP.pElement + ((size_t)IE) * (size_t)a_offset;
}

const double* ideal_observer::CDistribution::pBase_logP( int a_offset ) const
{
	if (a_offset<0 || a_offset>=IWH) return NULL;
	return ce_logP.pElement + ((size_t)IE) * (size_t)a_offset;
}

double* ideal_observer::CDistribution::pBase_realP( int a_offset )
{
	if (a_offset<0 || a_offset>=IWH) return NULL;
	return ce_realP.pElement + ((size_t)IE) * (size_t)a_offset;
}

const double* ideal_observer::CDistribution::pBase_realP( int a_offset ) const
{
	if (a_offset<0 || a_offset>=IWH) return NULL;
	return ce_realP.pElement + ((size_t)IE) * (size_t)a_offset;
}

void ideal_observer::CDistribution::Normalize_logP_and_Build_realP( int a_offset )
{
	int i, j, ia, ib;
	double* P_lp;
	double* P_rp;
	double ts;

	if (a_offset<0 || a_offset>=IWH) { ia = 0; ib = IWH-1; }
	else ia = ib = a_offset;
	for ( i = ia; i <= ib; ++i )
	{
		P_lp = pBase_logP(i);
		P_rp = pBase_realP(i);
		ts = -MAXIMUM_DOUBLE_NUM_V;
		for ( j = 0; j < IE; ++j )
			ts = CMyMath::log_Sum(ts, P_lp[j]);
		for ( j = 0; j < IE; ++j )
		{
			P_lp[j] = CMyMath::log_Product(P_lp[j], -ts);
			P_rp[j] = exp(P_lp[j]);
		}
	}
}

// [*] CObservation

ideal_observer::CObservation::CObservation()
{
	Fixation.Clear();
	Result.Clear();
}

ideal_observer::CObservation::CObservation(const CFixation& Fixation, const CSimpleTypeArray<int>& Result)
{
	this->Fixation = Fixation;
	this->Result = Result;
}

// [*] CObservationSet

ideal_observer::CObservationSet::CObservationSet()
{
	Clear();
}

void ideal_observer::CObservationSet::Clear()
{
	Size = 0;
	List.Clear();
}

void ideal_observer::CObservationSet::NewObservation(const CObservation& NewObservation)
{
	++Size;
	List.PushBack(NewObservation);
}

void ideal_observer::CObservationSet::_MakeFN_Observation( char* DstPtrChar, const char* strPN_ObservationBase, int ObservationIndex )
{
	sprintf(DstPtrChar, "%s%04d.txt", strPN_ObservationBase, ObservationIndex);
}

bool ideal_observer::CObservationSet::InputList( string strFN_SrcList, string strPN_SrcObservationBase, const CImageModelInstance& IMI )
{
	int i;
	char strTemp[10001];

	sprintf_s(strTemp, "Reading Observations from Disk...");
	CExternalTools::ConsoleTools::ClearCurrentLine();
	cout<<strTemp;

	Clear();
	ifstream inFile(strFN_SrcList.c_str());
	if (!(inFile>>Size) || Size<0)
	{
		Clear();
		cout<<" Failed;"<<endl;
		return false;
	}
	inFile.getline(strTemp, 10000);
	List.Resize(Size);
	for ( i = 0; i < Size; ++i )
	{
		CObservation& o = List.pElement[i];
		inFile.getline(strTemp, 10000);
		if (!(IMI.InputObservationFromText(o, strTemp)))
		{
			Clear();
			cout<<" Failed;"<<endl;
			return false;
		}
	}
	inFile.close();
	inFile.clear();
	cout<<" Done;"<<endl;
	return true;
}

void ideal_observer::CObservationSet::OutputList( string strFN_DstList, string strPN_DstObservationBase, const CImageModelInstance& IMI ) const
{
	int i;
	char strTemp[10001];

	//sprintf_s(strTemp, "Writing Observations to Disk...");
	//CExternalTools::ConsoleTools::ClearCurrentLine();
	//cout<<strTemp;

	CExternalTools::ConsoleTools::MakeDirectory(strPN_DstObservationBase, true);
	const char* ppp = strPN_DstObservationBase.c_str();

	ofstream outFile(strFN_DstList.c_str());
	outFile<<Size<<endl;
	for ( i = 0; i < Size; ++i )
	{
		const CObservation& o = List.pElement[i];
		_MakeFN_Observation(strTemp, ppp, i);
		outFile<<strTemp<<endl;
		IMI.OutputObservationToText(strTemp, o);
	}
	outFile.close();
	outFile.clear();

	//cout<<" Done;"<<endl;
}

// [*] CImageModelInstance::CAnalysisParameter

ideal_observer::CImageModelInstance::CAnalysisParameter::CAnalysisParameter()
{
	this->Sampling_IterationAmount = 131072;
	this->SeedAmount_Exponent = 0.5;
	this->SeedAmount_Coefficient = 1.0;
	this->T_AbandonThreshold = 0.001;
	this->Ra_AffectThreshold = 0.001;
	this->Cell_IndependentThreshold = 0.001;
}

// [*] CImageModelInstance

void ideal_observer::CImageModelInstance::InputObserveResultFromImage(
	CSimpleTypeArray<int>& matDstObserveResult,
	string strFN_SrcRGB,
	int IW, int IH, int CW, int CH)
{
	CImg<int> imgSrcRGB;
	CExternalTools::myLoadImage(strFN_SrcRGB, imgSrcRGB);
	InputObserveResultFromImage(matDstObserveResult, imgSrcRGB, IW, IH, CW, CH);
}

void ideal_observer::CImageModelInstance::InputObserveResultFromImage(
	CSimpleTypeArray<int>& matDstObserveResult,
	const CImg<int>& imgSrcRGB,
	int IW, int IH, int CW, int CH)
{
	matDstObserveResult.Clear();
	if (IW<1 || IH<1 || CW<3 || CH<3) return;

	CZoomableImage ZoomImage;
	CZoomConfig ZoomConfig;
	CImg<int> imgBlackWhite;
	int iiw, iih;
	int ix, iy, tr;
	CRectBox<int> iBox;
	CObjectEdglet toe;

	ZoomImage.myInit(imgSrcRGB);
	iiw = IW * CW;
	iih = IH * CH;
	ZoomConfig.myInit(
		1.0 * iiw / ZoomImage.W,
		1.0 * iih / ZoomImage.H);
	ZoomImage.DefaultColor.Fill(255.0);
	ZoomImage.BuildImage(imgBlackWhite, iiw, iih, ZoomConfig, CDenseImage::EXPANDING_MODE_COLOR);

	for ( iy = 0; iy < iih; ++iy ) for ( ix = 0; ix < iiw; ++ix )
		imgBlackWhite.at(ix, iy, 0) = (imgBlackWhite.at(ix, iy, 0) + imgBlackWhite.at(ix, iy, 1) + imgBlackWhite.at(ix, iy, 2)) / 3;

	//CImg<int> imgOutput = imgBlackWhite;
	//for ( iy = 0; iy < iih; ++iy ) for ( ix = 0; ix < iiw; ++ix )
	//{
	//	if (iy%CH!=0 && ix%CW!=0) continue;
	//	imgOutput.at(ix, iy, 0, 0) = 128;
	//	imgOutput.at(ix, iy, 0, 1) = 128;
	//	imgOutput.at(ix, iy, 0, 2) = 128;
	//}
	//CExternalTools::mySaveImage("ObserveResult.Show.jpg", imgOutput);
	////imgOutput.display();

	matDstObserveResult.Resize(IW * IH); matDstObserveResult.Clear();
	iBox.h = CH; iBox.w = CW;
	for ( iy = 0; iy < IH; ++iy ) for ( ix = 0; ix < IW; ++ix )
	{
		iBox.y = iy * CH;
		iBox.x = ix * CW;
		toe.y = iy;
		toe.x = ix;
		toe.edglet.FromImage(imgBlackWhite.data, iiw, iih, iBox);
		if (toe.edglet.fDoesExist())
			tr = toe.edglet.label;
		else
			tr = -1;
		matDstObserveResult.PushBack(tr);
	}
}

void ideal_observer::CImageModelInstance::OutputObserveResultToImage(
	CImg<int>& imgOutput,
	const CSimpleTypeArray<int>& matSrcObserveResult,
	int CW, int CH, int H_Extent, bool fDrawGrid) const
{
	if (IWH < 1) return;

	CImg<int> imgContext;
	CEdglet te;
	int i;
	CRectBox<int> iBox;
	int iiw, iih;
	iiw = IW * CW; iih = IH * CH;
	imgContext.assign(iiw, iih, 1, 1);
	imgContext.fill(255);
	iBox.w = CW; iBox.h = CH;
	for ( i = 0; i < IWH; ++i )
	{
		te.label = matSrcObserveResult.pElement[i];
		iBox.x = i%IW * CW;
		iBox.y = i/IW * CH;
		te.ToImage(imgContext.data, iiw, iih, iBox);
	}

	if (fDrawGrid)
		DrawGridImage(imgOutput, imgContext, IW, IH, CW, CH, H_Extent);
	else
	{
		imgOutput.assign(iiw, iih + H_Extent, 1, 3);
		imgOutput.fill(0);
		imgOutput.draw_image(imgContext, 0, 0, 0, 0);
		imgOutput.draw_image(imgContext, 0, 0, 0, 1);
		imgOutput.draw_image(imgContext, 0, 0, 0, 2);
	}
}

void ideal_observer::CImageModelInstance::OutputObserveResultToImage(
	CImg<int>& imgOutput,
	const CObservation& SrcObservation,
	int CW, int CH, int H_Extent, bool fDrawGrid) const
{
	OutputObserveResultToImage(imgOutput, SrcObservation.Result, CW, CH, H_Extent, fDrawGrid);
}

bool ideal_observer::CImageModelInstance::InputObservationFromText(
	CObservation& DstObservation,
	string strFN_SrcText) const
{
	DstObservation.Fixation.Clear();
	DstObservation.Result.Clear();
	if (IWH < 1) return false;

	ifstream inFile(strFN_SrcText.c_str());
	DstObservation.Result.Resize(IWH); DstObservation.Result.Fill(-1);
	int ix, iy, t;
	if (!(DstObservation.Fixation.InputFromStream(inFile))) return false;
	inFile>>ix>>iy; if (ix!=IW || iy!=IH) return false;
	for ( iy = 0; iy < IH; ++iy ) for ( ix = 0; ix < IW; ++ix )
	{
		if (!(inFile>>t) || (t<0) || (t>=IE)) continue; 
		DstObservation.Result.pElement[ix + iy * IW] = t;
	}
	inFile.close();
	inFile.clear();
	return true;
}

void ideal_observer::CImageModelInstance::OutputObservationToText(
	string strFN_DstText,
	const CObservation& SrcObservation) const
{
	if (IWH != (int)SrcObservation.Result.Size) return;

	ofstream outFile(strFN_DstText.c_str());
	int ix, iy, t;
	char strTemp[10001];
	SrcObservation.Fixation.OutputToStream(outFile);
	outFile<<"\t"<<IW<<"\t"<<IH<<endl;
	for ( iy = 0; iy < IH; ++iy )
	{
		for ( ix = 0; ix < IW; ++ix )
		{
			t = SrcObservation.Result.pElement[ix + iy * IW];
			sprintf_s(strTemp, " %02d", t);
			outFile<<strTemp;
		}
		outFile<<endl;
	}
	outFile.close();
	outFile.clear();
}

void ideal_observer::CImageModelInstance::InitObjectTemplate( string strFN_ObjectTemplate, int IW, int IH, int CW, int CH )
{
	this->IE = CEdglet::Amount_Edglet;

	this->IW = IW;
	this->IH = IH;
	this->IWH = IW * IH;
	InputObserveResultFromImage(this->matPrototypeD, strFN_ObjectTemplate, IW, IH, CW, CH);
	ObjectTemplate.myInit(this->matPrototypeD, IW, IH);

	this->_Transformation.Clear();
	this->_Fixation.Clear();

	this->InitModelParameters(0.75, 0.01, 0.99, 2.0);

	InitDefaultEdgletConfusionMatrix();
	InitDefaultEdgletPriorDistribution();
}

void ideal_observer::CImageModelInstance::InitModelParameters(double K1, double K2, double K3, double Sigma)
{
	K1 = (K1<0 ? 0 : (K1>1 ? 1 : K1));
	K2 = (K2<0 ? 0 : (K2>1 ? 1 : K2));
	K3 = (K3<0 ? 0 : (K3>1 ? 1 : K3));
	if (K2 > K3)
	{
		_ms = K2; K2 = K3; K3 = _ms;
	}
	this->K1 = K1;
	this->K2 = K2;
	this->K3 = K3;
	this->Sigma = Sigma;
	_ms = -0.5 / Sigma / Sigma;
}

void ideal_observer::CImageModelInstance::FixProbabilitySettings()
{
	int i, j;
	double* pd;
	double ss;

	if ((int)matr_Confusion.Size == IE * IE)
	{
		for ( i = 0; i < IE; ++i )
		{
			pd = matr_Confusion.pElement + i * IE;
			ss = 0.0;
			for ( j = 0; j < IE; ++j ) ss += pd[j];
			ss = (ss<EPS ? EPS : ss);
			ss = ss;
			for ( j = 0; j < IE; ++j ) pd[j] /= ss;
			ss = ss;
			ss = ss;
		}
	}

	if ((int)vect_EdgletPrior.Size == IE)
	{
		pd = vect_EdgletPrior.pElement;
		ss = 0.0;
		for ( j = 0; j < IE; ++j ) ss += pd[j];
		ss = (ss<EPS ? EPS : ss);
		ss = ss;
		for ( j = 0; j < IE; ++j ) pd[j] /= ss;
		ss = ss;
		ss = ss;
	}
}

void ideal_observer::CImageModelInstance::InitDefaultEdgletConfusionMatrix()
{
	matr_Confusion.Resize(IE * IE);
	matr_Confusion.Fill(1.0);
	CEdglet te;
	int i, ia, ib;
	int j, ja, jb;
	int tt;
	int t1, t2, t3, t4, t5;
	double* pd;
	tt = CEdglet::Amount_OriBin / 2;
	for ( i = 0; i < IE; ++i )
	{
		pd = matr_Confusion.pElement + i * IE;
		te.myInit(i); te.ToBinPair(ia, ib);
		for ( j = 0; j < IE; ++j )
		{
			te.myInit(j); te.ToBinPair(ja, jb);
			t1 = abs(ia-ja); if (t1 > tt) t1 = CEdglet::Amount_OriBin - t1;
			t2 = abs(ib-jb); if (t2 > tt) t2 = CEdglet::Amount_OriBin - t2;
			t3 = t1 + t2;
			t1 = abs(ib-ja); if (t1 > tt) t1 = CEdglet::Amount_OriBin - t1;
			t2 = abs(ia-jb); if (t2 > tt) t2 = CEdglet::Amount_OriBin - t2;
			t4 = t1 + t2;
			if (t3 < t4) t5 = t3; else t5 = t4;

			pd[j] = tt + tt - t5 + 1;
		}
	}
	FixProbabilitySettings();
}

void ideal_observer::CImageModelInstance::InitDefaultEdgletPriorDistribution()
{
	vect_EdgletPrior.Resize(IE);
	vect_EdgletPrior.Fill(1.0 / IE);
	FixProbabilitySettings();
}

void ideal_observer::CImageModelInstance::InitDefaultTransformationSet(CTransformationSet& DstTransformationSet)
{
	DstTransformationSet.Clear();
	CTransformation tT;
	int bx_min = 0 - ObjectTemplate.xLow;
	int bx_max = IW - ObjectTemplate.xHigh - 1;
	int by_min = 0 - ObjectTemplate.yLow;
	int by_max = IH - ObjectTemplate.yHigh - 1;
	int ix, iy;
	for ( iy = by_min; iy <= by_max; ++iy ) for ( ix = bx_min; ix <= bx_max; ++ix )
	{
		tT.myInit(ix, iy);
		DstTransformationSet.PushBack(tT, 1.0);
	}
	DstTransformationSet.NormalizeWeights();
}

int ideal_observer::CImageModelInstance::InitDefaultFixationSet(CFixationSet& DstFixationSet, double FixationStrideFactor)
{
	int ret = 1;

	DstFixationSet.Clear();
	CFixation tF;

	//int ix, iy;
	//for ( iy = 0; iy < IH; ++iy ) for ( ix = 0; ix < IW; ++ix )
	//{
	//	tF.myInit(ix, iy);
	//	DstFixationSet.PushBack(tF);
	//}

	double Radius_of_F = Sigma * FixationStrideFactor / sqrt(2.0);
	ret = CMyMath::TranslateNumber_Floor(Radius_of_F, (int)0);
	if (ret < 1) ret = 1;

	int ox = (IW-1) % ret / 2;
	int oy = (IH-1) % ret / 2;
	int ix, iy;
	for ( iy = oy; iy < IH; iy+=ret ) for ( ix = ox; ix < IW; ix+=ret )
	{
		tF.myInit(ix, iy);
		DstFixationSet.PushBack(tF);
	}

	return ret;
}

bool ideal_observer::CImageModelInstance::LoadEdgletConfusionMatrix(string strFN_EdgletConfusionMatrix)
{
	int i, j;
	double t;
	bool ret = true;
	ifstream inFile(strFN_EdgletConfusionMatrix.c_str());
	ret = ((inFile>>i) && (i==IE));
	matr_Confusion.Clear();
	for ( i = 0; ret && i < IE; ++i ) for ( j = 0; ret && j < IE; ++j )
	{
		ret = ret && (inFile>>t);
		t = (t<0 ? 0 : t);
		if (ret) matr_Confusion.PushBack(t);
	}
	inFile.close();
	inFile.clear();
	if (ret) FixProbabilitySettings();
	return ret;
}

void ideal_observer::CImageModelInstance::SaveEdgletConfusionMatrix(string strFN_EdgletConfusionMatrix)
{
	int i, j;
	char strTemp[10001];
	ofstream outFile(strFN_EdgletConfusionMatrix.c_str());
	outFile<<IE<<endl;
	for ( i = 0; i < IE; ++i )
	{
		for ( j = 0; j < IE; ++j )
		{
			sprintf_s(strTemp, "\t%10.8lf", matr_Confusion.pElement[j + i * IE]);
			outFile<<strTemp;
		}
		outFile<<endl;
	}
	outFile.close();
	outFile.clear();
}

bool ideal_observer::CImageModelInstance::LoadEdgletPriorDistribution(string strFN_EdgletPriorDistribution)
{
	int i;
	double t;
	bool ret = true;
	ifstream inFile(strFN_EdgletPriorDistribution.c_str());
	ret = ((inFile>>i) && (i==IE));
	vect_EdgletPrior.Clear();
	for ( i = 0; ret && i < IE; ++i )
	{
		ret = ret && (inFile>>t);
		t = (t<0 ? 0 : t);
		if (ret) vect_EdgletPrior.PushBack(t);
	}
	inFile.close();
	inFile.clear();
	if (ret) FixProbabilitySettings();
	return ret;
}

void ideal_observer::CImageModelInstance::SaveEdgletPriorDistribution(string strFN_EdgletPriorDistribution)
{
	int i;
	char strTemp[10001];
	ofstream outFile(strFN_EdgletPriorDistribution.c_str());
	outFile<<IE<<endl;
	for ( i = 0; i < IE; ++i )
	{
		sprintf_s(strTemp, "\t%10.8lf", vect_EdgletPrior.pElement[i]);
		outFile<<strTemp;
		outFile<<endl;
	}
	outFile.close();
	outFile.clear();
}

void ideal_observer::CImageModelInstance::SetTransformation(const CTransformation& T)
{
	this->_Transformation = T;
	this->_Refresh_PrototypeD_given_T();
}

void ideal_observer::CImageModelInstance::SetFixation(const CFixation& F)
{
	this->_Fixation = F;
}

const CTransformation& ideal_observer::CImageModelInstance::GetTransformation() const
{
	return _Transformation;
}

const CFixation& ideal_observer::CImageModelInstance::GetFixation() const
{
	return _Fixation;
}

void ideal_observer::CImageModelInstance::_Refresh_PrototypeD_given_T()
{
	ObjectTemplate.DoTransformation(matPrototypeD, _Transformation);
}

void ideal_observer::CImageModelInstance::Ouptut_PrototypeD_ToImage( CImg<int>& imgOutput, int CW, int CH, int H_Extent, bool fDrawGrid ) const
{
	OutputObserveResultToImage(imgOutput, matPrototypeD, CW, CH, H_Extent, fDrawGrid);
}

double ideal_observer::CImageModelInstance::func_z(int dx, int dy) const
{
	return K2 + (K3 - K2) * func_e(dx, dy);
}

double ideal_observer::CImageModelInstance::func_e(int dx, int dy) const
{
	return exp((dx*dx+dy*dy) * _ms);
}

double ideal_observer::CImageModelInstance::logP_Da_given_T( int a_offset, int d_a ) const
{
	double ret = -MAXIMUM_DOUBLE_NUM_V;

	if (a_offset<0 || a_offset>=IWH) return ret;
	//if (a_x<0 || a_x>=IW || a_y<0 || a_y>=IH) return ret;
	if (d_a<0 || d_a>=IE) return ret;

	int o_a = matPrototypeD.pElement[a_offset];

	if (o_a<0 || o_a>=IE)
		ret = vect_EdgletPrior.pElement[d_a];
	else
		ret = (d_a==o_a ? K1 : 0.0) + (1.0 - K1) * matr_Confusion.pElement[d_a + o_a * IE];
	ret = CMyMath::my_log(ret);
	return ret;
}

double ideal_observer::CImageModelInstance::logP_Ra_given_Da_F( int a_offset, int r_a, int d_a ) const
{
	double ret = -MAXIMUM_DOUBLE_NUM_V;

	if (a_offset<0 || a_offset>=IWH) return ret;
	//if (a_x<0 || a_x>=IW || a_y<0 || a_y>=IH) return ret;
	if (r_a<0 || r_a>=IE) return ret;
	if (d_a<0 || d_a>=IE) return ret;

	int dx = a_offset % IW - _Fixation.ob_x;
	int dy = a_offset / IW - _Fixation.ob_y;

	ret = func_z(dx, dy);
	ret = (r_a==d_a ? ret : 0.0) + (1.0 - ret) * matr_Confusion.pElement[r_a + d_a * IE];
	//ret = K2 * (1.0 - exp((dx*dx+dy*dy) * _ms));
	//ret = (r_a==d_a ? 1.0-ret : 0.0) + ret * matr_Confusion.pElement[r_a + d_a * IE];

	ret = CMyMath::my_log(ret);
	return ret;
}

void ideal_observer::CImageModelInstance::Predict_D_given_T( CDistribution& DstDistribution ) const
{
	DstDistribution.myInit(IW, IH);
	int itr_a, itr_d;
	double* P_d;
	for ( itr_a = 0; itr_a < IWH; ++itr_a )
	{
		P_d = DstDistribution.pBase_logP(itr_a);
		for ( itr_d = 0; itr_d < IE; ++itr_d )
			P_d[itr_d] = logP_Da_given_T(itr_a, itr_d);
	}
	DstDistribution.Normalize_logP_and_Build_realP();
}

void ideal_observer::CImageModelInstance::Predict_R_given_obsvD_F( CDistribution& DstDistribution, const CSimpleTypeArray<int>& obsv_D) const
{
	DstDistribution.myInit(IW, IH);
	int itr_a, itr_i, rsd;
	double* P_i;
	for ( itr_a = 0; itr_a < IWH; ++itr_a )
	{
		P_i = DstDistribution.pBase_logP(itr_a);
		rsd = obsv_D.pElement[itr_a];
		for ( itr_i = 0; itr_i < IE; ++itr_i )
			P_i[itr_i] = logP_Ra_given_Da_F(itr_a, itr_i, rsd);
	}
	DstDistribution.Normalize_logP_and_Build_realP();
}

void ideal_observer::CImageModelInstance::Predict_R_given_dtrbD_F(CDistribution& DstDistribution, const CDistribution& Src_dtrbD) const
{
	DstDistribution.myInit(IW, IH);
	int itr_a, itr_i, itr_d;
	double ts;
	double* P_i;
	const double* Q_d;
	for ( itr_a = 0; itr_a < IWH; ++itr_a )
	{
		P_i = DstDistribution.pBase_logP(itr_a);
		Q_d = Src_dtrbD.pBase_logP(itr_a);
		for ( itr_i = 0; itr_i < IE; ++itr_i )
		{
			P_i[itr_i] = -MAXIMUM_DOUBLE_NUM_V;
			for ( itr_d = 0; itr_d < IE; ++itr_d )
			{
				ts = logP_Ra_given_Da_F(itr_a, itr_i, itr_d);
				ts = CMyMath::log_Product(ts, Q_d[itr_d]);
				P_i[itr_i] = CMyMath::log_Sum(P_i[itr_i], ts);
			}
		}
	}
	DstDistribution.Normalize_logP_and_Build_realP();
}

void ideal_observer::CImageModelInstance::Predict_R_given_T_F(CDistribution& DstDistribution) const
{
	CDistribution dtrbD;
	this->Predict_D_given_T(dtrbD);
	this->Predict_R_given_dtrbD_F(DstDistribution, dtrbD);
}

double ideal_observer::CImageModelInstance::RandomSampleObserveResult(CSimpleTypeArray<int>& DstObserveResult, const CDistribution& SrcDistribution) const
{
	int i_a, j_R;
	double lsp, tlp;
	const double* Q_lp;
	const double* Q_rp;

	double ret = 0.0;
	DstObserveResult.Resize(IWH);
	DstObserveResult.Fill(-1);
	for ( i_a = 0; i_a < IWH; ++i_a )
	{
		lsp = -MAXIMUM_DOUBLE_NUM_V;
		Q_lp = SrcDistribution.pBase_logP(i_a);
		Q_rp = SrcDistribution.pBase_realP(i_a);

		j_R = CMyRandom::DiscreteSample(Q_rp, IE);
		tlp = Q_lp[j_R];
		ret = CMyMath::log_Product(ret, tlp);

		DstObserveResult.pElement[i_a] = j_R;
	}
	return ret;
}

double ideal_observer::CImageModelInstance::RandomSample_D_given_T(CSimpleTypeArray<int>& obsv_D) const
{
	CDistribution dtrbD;
	this->Predict_D_given_T(dtrbD);
	return RandomSampleObserveResult(obsv_D, dtrbD);
}

double ideal_observer::CImageModelInstance::RandomSample_R_given_obsvD_F(CObservation& DstObservation, const CSimpleTypeArray<int>& obsv_D) const
{
	CDistribution dtrbR;
	this->Predict_R_given_obsvD_F(dtrbR, obsv_D);
	DstObservation.Fixation = _Fixation;
	return RandomSampleObserveResult(DstObservation.Result, dtrbR);
}

double ideal_observer::CImageModelInstance::RandomSample_R_given_T_F(CObservation& DstObservation) const
{
	CDistribution dtrbR;
	this->Predict_R_given_T_F(dtrbR);
	DstObservation.Fixation = _Fixation;
	return RandomSampleObserveResult(DstObservation.Result, dtrbR);
}

double ideal_observer::CImageModelInstance::Analyse2_given_T_RR_FF(
	CDistribution& dtrb_D_given_T,
	CDistribution& dtrb_D_given_T_RR_FF,
	CDistribution& xyD_logP_RRa_given_Da_FF,
	CSimpleTypeArray<double>& mat_logP_RRa_given_T_FF,
	const CObservationSet& ObservationSet,
	const CTransformation& T)
{
	dtrb_D_given_T_RR_FF.myInit(IW, IH);
	xyD_logP_RRa_given_Da_FF.myInit(IW, IH);
	mat_logP_RRa_given_T_FF.Resize(IWH); mat_logP_RRa_given_T_FF.Fill(0.0);
	double logP_RR_given_T_FF = 0.0; // also return value

	double* P_logP_RRa_Da_given_T_FF;
	const double* Q_logP_Da_given_T;
	//double* Pd_logP_Da_given_T_RR_FF;
	double* PDa_logP_RRa_given_Da_FF;

	this->SetTransformation(T);
	this->Predict_D_given_T(dtrb_D_given_T);

	const CObservation* Q_Observation = ObservationSet.List.pElement;
	int i_a, i_F;
	int j_D;
	double tlp;

	logP_RR_given_T_FF = 0.0;
	for ( i_a = 0; i_a < IWH; ++i_a )
	{
		P_logP_RRa_Da_given_T_FF = dtrb_D_given_T_RR_FF.pBase_logP(i_a);
		Q_logP_Da_given_T = dtrb_D_given_T.pBase_logP(i_a);
		PDa_logP_RRa_given_Da_FF = xyD_logP_RRa_given_Da_FF.pBase_logP(i_a);
		mat_logP_RRa_given_T_FF.pElement[i_a] = -MAXIMUM_DOUBLE_NUM_V;
		for ( j_D = 0; j_D < IE; ++j_D )
		{
			double& tDa_logP_RRa_given_Da_FF = PDa_logP_RRa_given_Da_FF[j_D];
			tDa_logP_RRa_given_Da_FF = 0.0;
			for ( i_F = 0; i_F < ObservationSet.Size; ++i_F )
			{
				this->SetFixation(Q_Observation[i_F].Fixation);
				tlp = logP_Ra_given_Da_F(i_a, Q_Observation[i_F].Result.pElement[i_a], j_D);
				tDa_logP_RRa_given_Da_FF = CMyMath::log_Product(tDa_logP_RRa_given_Da_FF, tlp);
				if (!(tDa_logP_RRa_given_Da_FF > -MAXIMUM_DOUBLE_NUM_C)) break;
			}
			tlp = CMyMath::log_Product(tDa_logP_RRa_given_Da_FF, Q_logP_Da_given_T[j_D]);
			P_logP_RRa_Da_given_T_FF[j_D] = tlp;
			mat_logP_RRa_given_T_FF.pElement[i_a] = CMyMath::log_Sum(mat_logP_RRa_given_T_FF.pElement[i_a], tlp);
		}
		dtrb_D_given_T_RR_FF.Normalize_logP_and_Build_realP(i_a);
		logP_RR_given_T_FF = CMyMath::log_Product(logP_RR_given_T_FF, mat_logP_RRa_given_T_FF.pElement[i_a]);
	}

	return logP_RR_given_T_FF;
}

void ideal_observer::CImageModelInstance::Analyse2_given_TT_RR_FF_F(
	CSimpleTypeArray<double>& DstList_P_T, int& DstBest_T, double& DstH_T,
	CSimpleTypeArray<double>& DstList_CH_F, int& DstBest_F, CFixation& DstFixation,
	const CTransformationSet& SetT, int Answer_T_Index,
	const CObservationSet& ObservationSet,
	const CFixationSet& SetF,
	ostream& oLog,
	int Output_CW, int Output_CH, int Output_H_Extent,
	string strPN_Prediction_Output)
{
	//CMyRandom::InitRandomSeed(0);

	if (PRINT_ANALYSIS_IMAGE && strPN_Prediction_Output != "")
		CExternalTools::ConsoleTools::MakeDirectory(strPN_Prediction_Output, true);

	DstH_T = CMyMath::my_log(SetT.Size);
	DstList_P_T.Resize(SetT.Size);
	DstFixation.Clear();
	DstList_CH_F.Resize(SetF.Size);

	char strTemp[10001];
	CMyNormalizer_Sum normalizer_sum_log(true);
	int i_T, i_F, i_a, j_D, j_R, i_sample;
	double tlp, slp, ttlp, sslp;

	CSimpleTypeArray<double> vT_logP_T_Prior;
	CSimpleTypeArray<double> vT_realP_T_Prior;
	vT_logP_T_Prior.Resize(SetT.Size);
	vT_realP_T_Prior.Resize(SetT.Size);

	CSimpleTypeArray<double> vT_logP_T_Posterior;
	CSimpleTypeArray<double> vT_realP_T_Posterior;
	vT_logP_T_Posterior.Resize(SetT.Size);
	vT_realP_T_Posterior.Resize(SetT.Size);
	//CDistribution dtrb_D_Posterior;
	//dtrb_D_Posterior.myInit(IW, IH);

	// Setting Parameters
	int Amount_Sampling = AnalysisParameter.Sampling_IterationAmount;
	int Ideal_Amount_Seed = 1 + (int)(AnalysisParameter.SeedAmount_Coefficient * pow(1.0 * AnalysisParameter.Sampling_IterationAmount, AnalysisParameter.SeedAmount_Exponent));
	double AbandonThreshold_logP_T_Posterior = CMyMath::my_log(AnalysisParameter.T_AbandonThreshold / SetT.Size);
	double AffectThreshold_Ra = AnalysisParameter.Ra_AffectThreshold;
	double IndependentThreshold_Cell = AnalysisParameter.Cell_IndependentThreshold;

	sprintf_s(strTemp, "Amount of Sampling Iteration = %d; Ideal Amount of Seed = %d;",
		Amount_Sampling, Ideal_Amount_Seed);
	CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, true);
	CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);
	sprintf_s(strTemp, "Threshold : TA = %.4lf; RAM = %.4lf; IC = %.4lf;",
		AnalysisParameter.T_AbandonThreshold, AnalysisParameter.Ra_AffectThreshold, AnalysisParameter.Cell_IndependentThreshold);
	CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, true);
	CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);

	// Process Prior P(T)
	for ( i_T = 0; i_T < SetT.Size; ++i_T )	vT_logP_T_Prior.pElement[i_T] = CMyMath::my_log(SetT.Weights.pElement[i_T]);
	slp = normalizer_sum_log.Normalize(vT_logP_T_Prior);
	for ( i_T = 0; i_T < SetT.Size; ++i_T )	vT_realP_T_Prior.pElement[i_T] = exp(vT_logP_T_Prior.pElement[i_T]);

	CSimpleTypeArray<double> mat_aRD_logP_Ra_given_Da_F;
	mat_aRD_logP_Ra_given_Da_F.Resize(((size_t)IWH) * ((size_t)IE) * (size_t)IE);
	CSimpleTypeArray<CDistribution> vT_xyR_logP_RRa_Ra_given_T_FF_F;
	vT_xyR_logP_RRa_Ra_given_T_FF_F.Resize(SetT.Size);
	for ( i_T = 0; i_T < SetT.Size; ++i_T )
	{
		CDistribution& tT_xyR_logP_RRa_Ra_given_T_FF_F = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T];
		tT_xyR_logP_RRa_Ra_given_T_FF_F.myInit(IW, IH);
	}

	// Process Posterior Probabilities

	CSimpleTypeArray<CDistribution> vT_dtrb_D_given_T;
	vT_dtrb_D_given_T.Resize(SetT.Size);
	CSimpleTypeArray<CDistribution> vT_dtrb_D_given_T_RR_FF;
	vT_dtrb_D_given_T_RR_FF.Resize(SetT.Size);
	CDistribution xyD_logP_RRa_given_Da_FF;
	CSimpleTypeArray<CSimpleTypeArray<double>> vT_mat_logP_RRa_given_T_FF;
	vT_mat_logP_RRa_given_T_FF.Resize(SetT.Size);
	CSimpleTypeArray<double> vT_logP_RR_given_T_FF;
	vT_logP_RR_given_T_FF.Resize(SetT.Size);
	double logP_RR_given_FF;
	{
		logP_RR_given_FF = -MAXIMUM_DOUBLE_NUM_V;
		for ( i_T = 0; i_T < SetT.Size; ++i_T )
		{
			tlp = Analyse2_given_T_RR_FF(
				vT_dtrb_D_given_T.pElement[i_T],
				vT_dtrb_D_given_T_RR_FF.pElement[i_T],
				xyD_logP_RRa_given_Da_FF,
				vT_mat_logP_RRa_given_T_FF.pElement[i_T],
				ObservationSet, SetT.List.pElement[i_T]);
			vT_logP_RR_given_T_FF.pElement[i_T] = tlp;
			tlp = CMyMath::log_Product(
				tlp, vT_logP_T_Prior.pElement[i_T]);
			vT_logP_T_Posterior.pElement[i_T] = tlp;
			logP_RR_given_FF = CMyMath::log_Sum(logP_RR_given_FF, tlp);

			sprintf_s(strTemp, "Analyzed Posterior for each T [%d/%d]...", i_T, SetT.Size);
			CExternalTools::ConsoleTools::ClearCurrentLine();
			cout<<strTemp;
		} // i_T

		slp = normalizer_sum_log.Normalize(vT_logP_T_Posterior);
		DstH_T = 0.0;
		double Best_T_logP_Posterior = -MAXIMUM_DOUBLE_NUM_V;
		double Best_T_realP_Posterior = 0.0;
		for ( i_T = 0; i_T < SetT.Size; ++i_T )
		{
			tlp = exp(vT_logP_T_Posterior.pElement[i_T]);
			vT_realP_T_Posterior.pElement[i_T] = tlp;
			DstList_P_T.pElement[i_T] = tlp;
			DstH_T -= tlp * vT_logP_T_Posterior.pElement[i_T];

			if (Best_T_realP_Posterior < tlp)
			{
				Best_T_realP_Posterior = tlp;
				Best_T_logP_Posterior = vT_logP_T_Posterior.pElement[i_T];
				DstBest_T = i_T;
			}

			if (PRINT_ANALYSIS_IMAGE && strPN_Prediction_Output != "")
			{
				CImg<int> imgOutput;
				CImg<int> imgContext;

				const CDistribution& tT_dtrb_D_given_T = vT_dtrb_D_given_T.pElement[i_T];
				const CDistribution& tT_dtrb_D_given_T_RR_FF = vT_dtrb_D_given_T_RR_FF.pElement[i_T];
				//const CSimpleTypeArray<double>& tT_mat_logP_RRa_given_T_FF = vT_mat_logP_RRa_given_T_FF.pElement[i_T];
				//for ( i_a = 0; i_a < IWH; ++i_a )
				//{
				//	double* PD_logP_Da_given_T_RR_FF = dtrb_D_given_T_RR_FF.pBase_logP(i_a);
				//	const double* QD_logP_Da_given_T = tT_dtrb_D_given_T.pBase_logP(i_a);
				//	const double* QD_logP_RRa_given_Da_FF = xyD_logP_RRa_given_Da_FF.pBase_logP(i_a);
				//	for ( j_D = 0; j_D < IE; ++j_D )
				//	{
				//		tlp = QD_logP_Da_given_T[j_D];
				//		tlp = CMyMath::log_Product(tlp, QD_logP_RRa_given_Da_FF[j_D]);
				//		tlp = CMyMath::log_Product(tlp, -tT_mat_logP_RRa_given_T_FF.pElement[i_a]);
				//		PD_logP_Da_given_T_RR_FF[j_D] = tlp;
				//	}
				//	dtrb_D_given_T_RR_FF.Normalize_logP_and_Build_realP(i_a);
				//}

				tT_dtrb_D_given_T_RR_FF.OutputToImage(imgContext, Output_CW, Output_CH, 0, false);

				CSimpleTypeArray<CObjectEdglet> ListOE;
				this->ObjectTemplate.DoTransformation(ListOE, SetT.List.pElement[i_T], true);
				int i, tx, ty, ix, iy, tc;
				for ( i = 0; i < (int)ListOE.Size; ++i )
				{
					tx = ListOE.pElement[i].x;
					ty = ListOE.pElement[i].y;
					for ( iy = 0; iy < Output_CH; ++iy ) for ( ix = 0; ix < Output_CW; ++ix )
					{
						tc = imgContext.at(tx * Output_CW + ix, ty * Output_CH + iy, 0, 0);
						tc = (int)(tc * 0.800);
						imgContext.at(tx * Output_CW + ix, ty * Output_CH + iy, 0, 1) = tc;
						imgContext.at(tx * Output_CW + ix, ty * Output_CH + iy, 0, 2) = tc;
					}
				}
				DrawGridImage(imgOutput, imgContext, IW, IH, Output_CW, Output_CH, Output_H_Extent);
				sprintf_s(strTemp, "T = <%+04d,%+04d>; P(T|{I,F}) = %6.2lf%%;",
					SetT.List.pElement[i_T].base_x,
					SetT.List.pElement[i_T].base_y,
					DstList_P_T.pElement[i_T] * 100.0);
				CExternalTools::OutputTextToImage(imgOutput, strTemp, 4, imgOutput.height - Output_H_Extent + 2, 0, 0, 0, 255, 255, 255, 16);
				sprintf_s(strTemp, "%sD_given_T_R_F.%04d.jpg", strPN_Prediction_Output.c_str(), i_T);
				CExternalTools::mySaveImage(strTemp, imgOutput);
			} // (strPN_Prediction_Output != "")

			sprintf_s(strTemp, "T [%03d/%03d] = <%+04d,%+04d>, P = %6.2lf%% (%+16.8e)",
				i_T, SetT.Size,
				SetT.List.pElement[i_T].base_x,
				SetT.List.pElement[i_T].base_y,
				DstList_P_T.pElement[i_T] * 100.0,
				vT_logP_T_Posterior.pElement[i_T]);
			CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, false);
			oLog<<strTemp<<endl;
		} // i_T

		if (1) // Draw Distribution(T)
		{
			CImg<int> imgOutput;
			CImg<int> imgContext;
			int ix, iy, ciw, cih, ipx, ipy, tt;
			CImg<double> imgPT;
			imgPT.assign(IW, IH);
			imgPT.fill(0.0);
			for ( i_T = 0; i_T < SetT.Size; ++i_T )
			{
				ix = ObjectTemplate.xCenter + SetT.List.pElement[i_T].base_x;
				iy = ObjectTemplate.yCenter + SetT.List.pElement[i_T].base_y;
				if (ix<0 || ix>=IW || iy<0 || iy>=IH) continue;
				imgPT.at(ix, iy) = vT_realP_T_Posterior.pElement[i_T];
			}
			CExternalTools::NormalizeToGeneralIntGrayscaleImage(imgPT, imgOutput);
			ciw = Output_CW * IW;
			cih = Output_CH * IH;
			imgContext.assign(ciw, cih, 1, 3);
			for ( ipy = 0; ipy < cih; ++ipy ) for ( ipx = 0; ipx < ciw; ++ipx )
			{
				ix = ipx / Output_CW; iy = ipy / Output_CH;
				tt = imgOutput.at(ix, iy);
				imgContext.at(ipx, ipy, 0, 0) = tt;
				imgContext.at(ipx, ipy, 0, 1) = tt;
				imgContext.at(ipx, ipy, 0, 2) = tt;
			}
			DrawGridImage(imgOutput, imgContext, IW, IH, Output_CW, Output_CH, Output_H_Extent);
			sprintf_s(strTemp, "H(T) = %.4lf; bT = <%+04d,%+04d>; P(bT|{I,F}) = %6.2lf%%;",
				DstH_T,
				SetT.List.pElement[DstBest_T].base_x,
				SetT.List.pElement[DstBest_T].base_y,
				DstList_P_T.pElement[DstBest_T] * 100.0);
			CExternalTools::OutputTextToImage(imgOutput, strTemp, 4, imgOutput.height - Output_H_Extent + 2, 0, 0, 0, 255, 255, 255, 16);
			sprintf_s(strTemp, "%sT_given_R_F.jpg", strPN_Prediction_Output.c_str());
			CExternalTools::mySaveImage(strTemp, imgOutput);
		} // Draw Distribution(T)

		sprintf_s(strTemp, "Current H(T) = %+8.4lf (%+16.8e);", DstH_T, log(DstH_T));
		CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, true);
		oLog<<strTemp<<endl;
		sprintf_s(strTemp, "Best T [%03d/%03d] = <%+04d,%+04d>, P = %6.2lf%% (%+16.8e);",
			DstBest_T, SetT.Size,
			SetT.List.pElement[DstBest_T].base_x,
			SetT.List.pElement[DstBest_T].base_y,
			DstList_P_T.pElement[DstBest_T] * 100.0,
			vT_logP_T_Posterior.pElement[DstBest_T]);
		CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, true);
		oLog<<strTemp<<endl;
	} // Process Posterior Probabilities

	// Filtering Unnecessary T
	CSimpleTypeArray<int> Waiting_T_List;
	int i_i_T;
	{
		Waiting_T_List.Resize(SetT.Size); Waiting_T_List.Clear();
		for ( i_T = 0; i_T < SetT.Size; ++i_T )
		{
			if (vT_logP_T_Posterior.pElement[i_T] < AbandonThreshold_logP_T_Posterior)
			{
				continue;
			}
			Waiting_T_List.PushBack(i_T);
		}
		sprintf_s(strTemp, "Selected |T| = %d;", (int)Waiting_T_List.Size);
		CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, true);
		CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);
	}

	// Process CH[F] for each F

	int FL, FH;
	FL = FH = -1;
	double CH_F_L, CH_F_H;
	double CH_1, CH_2;
	string stc_1, stc_2, stc_F;
	CTimer timCH_1, timCH_2, timCH_F, timCH_sfs;
	timCH_1.Reset(); timCH_2.Reset(); timCH_F.Reset(); timCH_sfs.Reset();

	sprintf_s(strTemp, "F_[i/amount]\tF_<fx,fy>\tCH_F.First\tCH_F.Second\tCH_F.Full\tbF_[i]\tbF_<fx,fy>\tCH_bF.Full");
	CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);

	for ( i_F = 0; i_F < SetF.Size; ++i_F )
	{ 
		double& CH_F = DstList_CH_F.pElement[i_F];
		CH_1 = CH_2 = 0.0;

		timCH_1.Start(); timCH_F.Start();

		if (1) // Process P(RR,R|T,FF,F)
		{
			double* PR_logP_RRa_Ra_given_T_FF_F;
			
			this->SetFixation(SetF.List.pElement[i_F]);
			for ( i_a = 0; i_a < IWH; ++i_a )
			{
				double* _Pa = mat_aRD_logP_Ra_given_Da_F.pElement + ((size_t)i_a) * ((size_t)IE) * (size_t)IE;
				for ( j_R = 0; j_R < IE; ++j_R )
				{
					double* _PaR = _Pa + ((size_t)j_R) * (size_t)IE;
					for ( j_D = 0; j_D < IE; ++j_D )
						_PaR[j_D] = this->logP_Ra_given_Da_F(i_a, j_R, j_D);
				} // j_R
			} // i_a

			for ( i_T = 0; i_T < SetT.Size; ++i_T )
			{
				CDistribution& tT_xyR_logP_RRa_Ra_given_T_FF_F = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T];
				const CDistribution& tT_dtrb_D_given_T = vT_dtrb_D_given_T.pElement[i_T];

				//double ss_max_p_p, ss_max_a_p;
				//ss_max_p_p = 0.0;
				//double ss_min_p_p, ss_min_a_p;
				//ss_min_p_p = 0.0;

				for ( i_a = 0; i_a < IWH; ++i_a )
				{
					PR_logP_RRa_Ra_given_T_FF_F = tT_xyR_logP_RRa_Ra_given_T_FF_F.pBase_logP(i_a);
					const double* QD_logP_RRa_given_Da_FF = xyD_logP_RRa_given_Da_FF.pBase_logP(i_a);
					const double* QD_D_given_T = tT_dtrb_D_given_T.pBase_logP(i_a);
					const double* _Qa = mat_aRD_logP_Ra_given_Da_F.pElement + ((size_t)i_a) * ((size_t)IE) * (size_t)IE;

					//double avg_p, avg_pp, var_pp, rel_var_p;
					//avg_p = -MAXIMUM_DOUBLE_NUM_V;
					//avg_pp = -MAXIMUM_DOUBLE_NUM_V;

					for ( j_R = 0; j_R < IE; ++j_R )
					{
						double& logP_RRa_Ra_given_T_FF_F = PR_logP_RRa_Ra_given_T_FF_F[j_R];
						logP_RRa_Ra_given_T_FF_F = -MAXIMUM_DOUBLE_NUM_V;
						const double* _QaR = _Qa + ((size_t)j_R) * (size_t)IE;
						for ( j_D = 0; j_D < IE; ++j_D )
						{
							tlp = _QaR[j_D]; // P(Ra|Da,F)
							tlp = CMyMath::log_Product(tlp, QD_logP_RRa_given_Da_FF[j_D]); // P(RRa|Da,FF)
							tlp = CMyMath::log_Product(tlp, QD_D_given_T[j_D]); // P(Da|T)
							logP_RRa_Ra_given_T_FF_F = CMyMath::log_Sum(logP_RRa_Ra_given_T_FF_F, tlp);
						} // j_D

						//tlp = logP_RRa_Ra_given_T_FF_F;
						//avg_p = CMyMath::log_Sum(avg_p, tlp);
						//tlp = CMyMath::log_Product(tlp, tlp);
						//avg_pp = CMyMath::log_Sum(avg_pp, tlp);

						//if (j_R==0 || ss_max_a_p<logP_RRa_Ra_given_T_FF_F)
						//	ss_max_a_p = logP_RRa_Ra_given_T_FF_F;
						//if (j_R==0 || ss_min_a_p>logP_RRa_Ra_given_T_FF_F)
						//	ss_min_a_p = logP_RRa_Ra_given_T_FF_F;
					} // j_R

					//avg_p = CMyMath::log_Product(avg_p, -log(IE * 1.0));
					//avg_pp = CMyMath::log_Product(avg_pp, -log(IE * 1.0));
					//tlp = CMyMath::log_Product(avg_p, avg_p);
					//var_pp = CMyMath::log_Minus(avg_pp, tlp);
					//tlp = CMyMath::log_Power(var_pp, 0.5);	
					//rel_var_p = CMyMath::log_Product(tlp, -avg_p);
					//tlp = exp(rel_var_p);

					//ss_max_a_p = CMyMath::log_Product(ss_max_a_p, -vT_mat_logP_RRa_given_T_FF.pElement[i_T].pElement[i_a]);
					//ss_max_p_p = CMyMath::log_Product(ss_max_p_p, ss_max_a_p);
					//ss_min_a_p = CMyMath::log_Product(ss_min_a_p, -vT_mat_logP_RRa_given_T_FF.pElement[i_T].pElement[i_a]);
					//ss_min_p_p = CMyMath::log_Product(ss_min_p_p, ss_min_a_p);
				} // i_a

				//sprintf_s(strTemp, "%10.8lf  %10.8lf  %10.8lf",
				//	vT_realP_T_Posterior.pElement[i_T], ss_min_p_p, ss_max_p_p);
				//CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);
			} // i_T
		} // Process P(RR,R|T,FF,F)

		if (1) // Make CH_1
		{
			sprintf_s(strTemp, "F [%d/%d]; Computing CH_1...", i_F, (int)SetF.Size);
			CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, false);
			CH_1 = 0.0;
			for ( i_T = 0; i_T < SetT.Size; ++i_T )
			{
				const CDistribution& tT_xyR_logP_RRa_Ra_given_T_FF_F = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T];
				const CSimpleTypeArray<double>& tT_mat_logP_RRa_given_T_FF = vT_mat_logP_RRa_given_T_FF.pElement[i_T];
				if (vT_logP_T_Posterior.pElement[i_T] < -MAXIMUM_DOUBLE_NUM_C) continue;
				sslp = 0.0;
				for ( i_a = 0; i_a < IWH; ++i_a )
				{
					slp = 0.0;
					const double* QR_logP_RRa_Ra_given_T_FF_F = tT_xyR_logP_RRa_Ra_given_T_FF_F.pBase_logP(i_a);
					double logP_RRa_given_T_FF = tT_mat_logP_RRa_given_T_FF.pElement[i_a];
					if (logP_RRa_given_T_FF < -MAXIMUM_DOUBLE_NUM_C) continue;
					for ( j_R = 0; j_R < IE; ++j_R )
					{
						tlp = QR_logP_RRa_Ra_given_T_FF_F[j_R];
						if (tlp < -MAXIMUM_DOUBLE_NUM_C) continue;
						ttlp = CMyMath::log_Product(tlp, -logP_RRa_given_T_FF);
						slp += exp(ttlp) * tlp;
					}
					sslp += slp;
				}
				sslp += vT_logP_T_Prior.pElement[i_T];
				ttlp = vT_realP_T_Posterior.pElement[i_T];
 				CH_1 -= ttlp * sslp;
			}
		} // Make CH_1

		timCH_1.Stop(); timCH_2.Start();

#if 0
		if (1) // Make CH_2 - Alan
		{
			double to_divide;
			int sampled_T;
			CSimpleTypeArray<int> obsv_D;
			obsv_D.Resize(IWH);
			CObservation oR_given_D_F;
			oR_given_D_F.Fixation = SetF.List.pElement[i_F];
			oR_given_D_F.Result.Resize(IWH);

			timCH_sfs.Reset(); timCH_sfs.Start();

			CSimpleTypeArray<double> wList;
			CSimpleTypeArray<double> vList;
			wList.Clear(); vList.Clear();

			CH_2 = -MAXIMUM_DOUBLE_NUM_V;
			to_divide = -MAXIMUM_DOUBLE_NUM_V;
			for ( sampled_T = 0; sampled_T < SetT.Size; ++sampled_T )
			{
				double w_w = 0.0;
				double w_B = 0.0;
				double w_F = 0.0;

				w_w = vT_logP_T_Prior.pElement[sampled_T];
				//w_w = CMyMath::log_Product(w_w, -logP_RR_given_FF);

				this->SetTransformation(SetT.List.pElement[sampled_T]);
				CSimpleTypeArray<int> tOT;
				this->ObjectTemplate.DoTransformation(tOT, this->_Transformation);
				CSimpleTypeArray<int> z_F; z_F.Clear();
				CSimpleTypeArray<int> z_B; z_B.Clear();
				for ( i_a = 0; i_a < IWH; ++i_a )
				{
					if (tOT.pElement[i_a] < 0)
						z_B.PushBack(i_a);
					else
						z_F.PushBack(i_a);
				}

				for ( i_sample = 0; i_sample < Amount_Sampling; ++i_sample )
				{
					const CDistribution& tT_dtrb_D_given_T = vT_dtrb_D_given_T.pElement[sampled_T];
					size_t offset_a = ((size_t)IE) * (size_t)IE;
					size_t offset_Ra = (size_t)IE;
					CSimpleTypeArray<double> tv_realP;
					tv_realP.Resize(IE);

					CSimpleTypeArray<double> vT_logP_RR_R_T_given_FF_F;
					vT_logP_RR_R_T_given_FF_F = vT_logP_T_Prior;

					CSimpleTypeArray<double> vT_FG_logP_RR_R_T_given_FF_F;
					vT_FG_logP_RR_R_T_given_FF_F.Resize(SetT.Size);
					vT_FG_logP_RR_R_T_given_FF_F.Fill(0.0);
					CSimpleTypeArray<double> vT_BG_logP_RR_R_T_given_FF_F;
					vT_BG_logP_RR_R_T_given_FF_F.Resize(SetT.Size);
					vT_BG_logP_RR_R_T_given_FF_F.Fill(0.0);

					int ii_a_B, ii_a_F;

					w_B = 0.0;
					for ( ii_a_B = 0; ii_a_B < (int)z_B.Size; ++ii_a_B )
					{
						i_a = z_B.pElement[ii_a_B];

						int& t_Da = obsv_D.pElement[i_a];
						int& t_Ra = oR_given_D_F.Result.pElement[i_a];

						t_Da = CMyRandom::DiscreteSample(tT_dtrb_D_given_T.pBase_realP(i_a), IE);

						const double* QR_logP_Ra_given_Da_F = mat_aRD_logP_Ra_given_Da_F.pElement + offset_a * (size_t)i_a + (size_t)t_Da;
						for ( j_R = 0; j_R < IE; ++j_R ) tv_realP.pElement[j_R] = exp(QR_logP_Ra_given_Da_F[offset_Ra * (size_t)j_R]);
						t_Ra = CMyRandom::DiscreteSample(tv_realP.pElement, IE);

						tlp = xyD_logP_RRa_given_Da_FF.pBase_logP(i_a)[t_Da];
						w_B = CMyMath::log_Product(w_B, tlp);

						for ( i_T = 0; i_T < SetT.Size; ++i_T )
						{
							double& logP_RR_R_T_given_FF_F = vT_logP_RR_R_T_given_FF_F.pElement[i_T];
							ttlp = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T].pBase_logP(i_a)[t_Ra];
							logP_RR_R_T_given_FF_F = CMyMath::log_Product(logP_RR_R_T_given_FF_F, ttlp);
							vT_BG_logP_RR_R_T_given_FF_F.pElement[i_T] = CMyMath::log_Product(vT_BG_logP_RR_R_T_given_FF_F.pElement[i_T], ttlp);
						}
					}

					w_F = 0.0;
					for ( ii_a_F = 0; ii_a_F < (int)z_F.Size; ++ii_a_F )
					{
						i_a = z_F.pElement[ii_a_F];

						int& t_Da = obsv_D.pElement[i_a];
						int& t_Ra = oR_given_D_F.Result.pElement[i_a];

						const double* QQ = tT_dtrb_D_given_T.pBase_realP(i_a);
						t_Da = CMyRandom::DiscreteSample(QQ, IE);

						const double* QR_logP_Ra_given_Da_F = mat_aRD_logP_Ra_given_Da_F.pElement + offset_a * (size_t)i_a + (size_t)t_Da;
						for ( j_R = 0; j_R < IE; ++j_R ) tv_realP.pElement[j_R] = exp(QR_logP_Ra_given_Da_F[offset_Ra * (size_t)j_R]);
						t_Ra = CMyRandom::DiscreteSample(tv_realP.pElement, IE);

						tlp = xyD_logP_RRa_given_Da_FF.pBase_logP(i_a)[t_Da];
						w_F = CMyMath::log_Product(w_F, tlp);

						for ( i_T = 0; i_T < SetT.Size; ++i_T )
						{
							double& logP_RR_R_T_given_FF_F = vT_logP_RR_R_T_given_FF_F.pElement[i_T];
							ttlp = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T].pBase_logP(i_a)[t_Ra];
							logP_RR_R_T_given_FF_F = CMyMath::log_Product(logP_RR_R_T_given_FF_F, ttlp);
							vT_FG_logP_RR_R_T_given_FF_F.pElement[i_T] = CMyMath::log_Product(vT_FG_logP_RR_R_T_given_FF_F.pElement[i_T], ttlp);
						}
					}

					sslp = -MAXIMUM_DOUBLE_NUM_V;
					for ( i_T = 0; i_T < SetT.Size; ++i_T )
						sslp = CMyMath::log_Sum(sslp, vT_logP_RR_R_T_given_FF_F.pElement[i_T]);

					slp = w_w;
					slp = CMyMath::log_Product(slp, w_B);
					slp = CMyMath::log_Product(slp, w_F);
					if (slp > -MAXIMUM_DOUBLE_NUM_C)
						slp = slp;

					tlp = CMyMath::my_log(-sslp);
					tlp = CMyMath::log_Product(tlp, slp);
					CH_2 = CMyMath::log_Sum(CH_2, tlp);
					to_divide = CMyMath::log_Sum(to_divide, slp);

					wList.PushBack(slp);
					vList.PushBack(tlp);
					sprintf_s(strTemp, "%06d\t\t%+12.6lf\t%+12.6lf\t%+12.6lf\t\t%+12.6lf\t%+12.6lf\t\t%+12.6lf\t%+12.6lf",
						i_sample+1,
						w_w, w_B, w_F,
						vT_BG_logP_RR_R_T_given_FF_F.pElement[0], vT_FG_logP_RR_R_T_given_FF_F.pElement[0], 
						slp, sslp);
					if (CMyRandom::NextInt(100) == 0)
						CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);

					if (CMyRandom::NextInt(Amount_Sampling/100) == 0)
					{
						timCH_sfs.Stop();
						stc_2 = CExternalTools::GetTimeText(
							timCH_sfs.GetCumulativeTime() * (Amount_Sampling - i_sample - 1) / (i_sample + 1));
						timCH_sfs.Start();

						tlp = CMyMath::log_Product(CH_2, -to_divide);
						tlp = -exp(tlp);
						sprintf_s(strTemp, "Analyzing F [%d/%d]; Progress [%02.0lf%%], CH = %+8.4lf; Time Left = %s;",
							i_F, SetF.Size, 100.0*(i_sample+1)/Amount_Sampling, CH_1 + tlp, stc_2.c_str());
						CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, false);

						//sprintf_s(strTemp, "Analyzing F [%d/%d]; Progress [%02.0lf%%], Times =\t%d\t, CH =\t%+9.6lf\t; Time Left = %s;",
						//	i_F, SetF.Size, 100.0*(i_sample+1)/Amount_Sampling, (i_sample+1), CH_1 + tlp, stc_2.c_str());
						//CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);
					}
				} // i_sample
			} // sampled_T

			timCH_sfs.Stop();

			double SCALE_STEP = 1.2;
			int BASE_NUMBER = 100;
			int i_scale, i_times;
			int N_Data = (int)vList.Size;
			for ( i_scale = 0; ; ++i_scale )
			{
				int ns = (int)(BASE_NUMBER * pow(SCALE_STEP, i_scale));
				if (ns > N_Data)
					ns = N_Data;
				
				for ( i_times = 0; i_times < 100; ++i_times )
				{
					int iA = ns;
					int iB = N_Data;
					sslp = -MAXIMUM_DOUBLE_NUM_V;
					slp = -MAXIMUM_DOUBLE_NUM_V;
					for ( i_sample = 0; i_sample < N_Data; ++i_sample, --iB )
					{
						if (CMyRandom::NextInt(iB) >= iA) continue;

						--iA;
						sslp = CMyMath::log_Sum(sslp, vList.pElement[i_sample]);
						slp = CMyMath::log_Sum(slp, wList.pElement[i_sample]);
					}
					sslp = CMyMath::log_Product(sslp, -slp);
					sslp = -exp(sslp);
					sprintf_s(strTemp, "%6d\t%12.6lf", ns, sslp);
					CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);
				}
				if (ns >= N_Data * 0.60) break;
			}

			CH_2 = CMyMath::log_Product(CH_2, -to_divide);
			CH_2 = -exp(CH_2);
		}
#else

		CSimpleTypeArray<int> Waiting_Cell_List;
		double Ideal_T_Count;
		double Selected_T_Count;

		if (1) // Make CH_2
		{
			double to_divide;

			sprintf_s(strTemp, "F [%d/%d]; Computing CH_2...", i_F, (int)SetF.Size);
			CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, false);

			CSimpleTypeArray<int> obsv_D;
			obsv_D.Resize(IWH);
			CObservation oR_given_D_F;
			oR_given_D_F.Fixation = SetF.List.pElement[i_F];
			oR_given_D_F.Result.Resize(IWH);
			int sampled_T;

			// Making P(R|T,RR,FF,F)
			// Also Decide Independent Cells
			CSimpleTypeArray<CDistribution> vT_xyR_logP_Ra_given_T_RRa_FF_F;
			int i_i_a;
			CSimpleTypeArray<double> Extra_log_ppT_List;
			//double Extra_CH_2;
			{
				CDistribution MaxValue_xyR_logP_Ra_given_T_RRa_FF_F;
				MaxValue_xyR_logP_Ra_given_T_RRa_FF_F.myInit(IW, IH);
				CDistribution MinValue_xyR_logP_Ra_given_T_RRa_FF_F;
				MinValue_xyR_logP_Ra_given_T_RRa_FF_F.myInit(IW, IH);

				vT_xyR_logP_Ra_given_T_RRa_FF_F = vT_xyR_logP_RRa_Ra_given_T_FF_F;
				for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
				{
					i_T = Waiting_T_List.pElement[i_i_T];
				//for ( i_T = 0; i_T < SetT.Size; ++i_T )
				//{
					CDistribution& tT_xyR_logP_Ra_given_T_RRa_FF_F = vT_xyR_logP_Ra_given_T_RRa_FF_F.pElement[i_T];
					tT_xyR_logP_Ra_given_T_RRa_FF_F.Normalize_logP_and_Build_realP();
					for ( i_a = 0; i_a < IWH; ++i_a )
					{
						const double* Q_T_Value = tT_xyR_logP_Ra_given_T_RRa_FF_F.pBase_logP(i_a);
						double* PMaxValue = MaxValue_xyR_logP_Ra_given_T_RRa_FF_F.pBase_logP(i_a);
						double* PMinValue = MinValue_xyR_logP_Ra_given_T_RRa_FF_F.pBase_logP(i_a);
						for ( j_R = 0; j_R < IE; ++j_R )
						{
							if (i_T==0 || PMaxValue[j_R]<Q_T_Value[j_R]) PMaxValue[j_R] = Q_T_Value[j_R];
							if (i_T==0 || PMinValue[j_R]>Q_T_Value[j_R]) PMinValue[j_R] = Q_T_Value[j_R];
						}
					}
				}

				//Extra_CH_2 = 0.0;
				Waiting_Cell_List.Clear();
				Extra_log_ppT_List = vT_logP_T_Prior;
				for ( i_a = 0; i_a < IWH; ++i_a )
				{
					const double* QMaxValue = MaxValue_xyR_logP_Ra_given_T_RRa_FF_F.pBase_logP(i_a);
					const double* QMinValue = MinValue_xyR_logP_Ra_given_T_RRa_FF_F.pBase_logP(i_a);
					tlp = 0.0;
					for ( j_R = 0; j_R < IE; ++j_R )
					{
						ttlp = QMaxValue[j_R] - QMinValue[j_R];
						if (tlp < ttlp)
							tlp = ttlp;
					}
					//tlp = 1.00;
					if (tlp > IndependentThreshold_Cell)
						Waiting_Cell_List.PushBack(i_a);
					else
					{
						i_i_T = 0; i_T = Waiting_T_List.pElement[i_i_T];
						// i_T = 0;
						{
							sslp = 0.0;
							const double* Q_log_T_Value = vT_xyR_logP_Ra_given_T_RRa_FF_F.pElement[i_T].pBase_logP(i_a);
							const double* Q_real_T_Value = vT_xyR_logP_Ra_given_T_RRa_FF_F.pElement[i_T].pBase_realP(i_a);
							for ( j_R = 0; j_R < IE; ++j_R )
							{
								sslp += Q_real_T_Value[j_R] * Q_log_T_Value[j_R];
							}
							//Extra_CH_2 = CMyMath::log_Product(Extra_CH_2, sslp);
						}

						//j_R = 0;
						{
							//ttlp = vT_xyR_logP_Ra_given_T_RRa_FF_F.pElement[i_T].pBase_logP(i_a)[j_R];
							for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
							{
								i_T = Waiting_T_List.pElement[i_i_T];
							//for ( i_T = 0; i_T < SetT.Size; ++i_T )
							//{
								double& tT_ppT = Extra_log_ppT_List.pElement[i_T];
								
								tlp = vT_mat_logP_RRa_given_T_FF.pElement[i_T].pElement[i_a];
								tlp = CMyMath::log_Product(tlp, sslp);
								tT_ppT = CMyMath::log_Product(tT_ppT, tlp);
							}
						}
					}
				}
			}

			// Making Affecting Matrix;
			size_t AM_Size_IE_IE = ((size_t)IE) * (size_t)IE;
			size_t AM_Size_IWH_IE_IE = AM_Size_IE_IE * (size_t)IWH;
			CSimpleTypeArray<int> AffectMatrix_Amount_T;
			CSimpleTypeArray<CSimpleTypeArray<int>> AffectMatrix_T;
			CSimpleTypeArray<CSimpleTypeArray<double>> AffectMatrix_Delta_logP;
			CSimpleTypeArray<CSimpleTypeArray<int>> AffectMatrix_Delta_ZeroCount;
			{
				//sprintf_s(strTemp, "F [%d/%d]; Making Affecting Matrix...", i_F, (int)SetF.Size);
				//CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, false);

				AffectMatrix_Amount_T.Resize(AM_Size_IWH_IE_IE);
				AffectMatrix_Amount_T.Fill(-1);
				AffectMatrix_T.Resize(AM_Size_IWH_IE_IE);
				AffectMatrix_Delta_logP.Resize(AM_Size_IWH_IE_IE);
				AffectMatrix_Delta_ZeroCount.Resize(AM_Size_IWH_IE_IE);
				//int old_Ra, t_Ra;
				//for ( i_a = 0; i_a < IWH; ++i_a ) for ( old_Ra = 0; old_Ra < IE; ++old_Ra ) for ( t_Ra = 0; t_Ra < IE; ++t_Ra )
				//{
				//	int& ta_Amount_T = AffectMatrix_Amount_T.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a];
				//	CSimpleTypeArray<int>& ta_List_T = AffectMatrix_T.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a];
				//	CSimpleTypeArray<double>& ta_List_Delta_logP = AffectMatrix_Delta_logP.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a];
				//	CSimpleTypeArray<int>& ta_List_Delta_ZeroCount = AffectMatrix_Delta_ZeroCount.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a];
				//	ta_Amount_T = 0; ta_List_T.Clear(); ta_List_Delta_logP.Clear(); ta_List_Delta_ZeroCount.Clear();

				//	double t_logP;
				//	int t_ZeroCount;

				//	for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
				//	{
				//		i_T = Waiting_T_List.pElement[i_i_T];
				//	//for ( i_T = 0; i_T < SetT.Size; ++i_T )
				//	//{
				//		const double* QRa_xyR_logP_RRa_Ra_given_T_FF_F = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T].pBase_logP(i_a);
				//		t_logP = 0.0; t_ZeroCount = 0;

				//		ttlp = QRa_xyR_logP_RRa_Ra_given_T_FF_F[old_Ra];
				//		if (ttlp > -MAXIMUM_DOUBLE_NUM_C)
				//			t_logP = CMyMath::log_Product(t_logP, -ttlp);
				//		else
				//			--t_ZeroCount;

				//		ttlp = QRa_xyR_logP_RRa_Ra_given_T_FF_F[t_Ra];
				//		if (ttlp > -MAXIMUM_DOUBLE_NUM_C)
				//			t_logP = CMyMath::log_Product(t_logP, +ttlp);
				//		else
				//			++t_ZeroCount;

				//		if (t_ZeroCount==0 && abs(t_logP)<AffectThreshold_Ra)
				//			continue;

				//		if (abs(t_logP) < 0.01)
				//			t_logP = t_logP;

				//		++ta_Amount_T;
				//		ta_List_T.PushBack(i_T);
				//		ta_List_Delta_logP.PushBack(t_logP);
				//		ta_List_Delta_ZeroCount.PushBack(t_ZeroCount);
				//	}
				//}
			}

			//for ( int i_try = 0; i_try < 20; ++i_try )
			//if (1) // Gibbs Sampling
			{
				CH_2 = 0.0;
				to_divide = 0.0;

				CSimpleTypeArray<double> vT_logP_RR_R_T_given_FF_F;
				vT_logP_RR_R_T_given_FF_F.Resize(SetT.Size);
				CSimpleTypeArray<int> vT_CountZeroP_RR_R_T_given_FF_F;
				vT_CountZeroP_RR_R_T_given_FF_F.Resize(SetT.Size);
				//CSimpleTypeArray<double> vT_realP_T_given_RR_R_FF_F;
				//vT_realP_T_given_RR_R_FF_F.Resize(SetT.Size); vT_realP_T_given_RR_R_FF_F.Fill(0.0);

				timCH_sfs.Reset(); timCH_sfs.Start();

				//CSimpleTypeArray<double> vList;
				//vList.Clear();

				// Preparing Seeds
				CSimpleTypeArray<int> Ideal_T_List;
				CSimpleTypeArray<bool> fMakeNewSeed_List;
				{
					Ideal_T_List.Resize(Amount_Sampling);
					fMakeNewSeed_List.Resize(SetT.Size); fMakeNewSeed_List.Fill(false);
					for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
					{
						i_T = Waiting_T_List.pElement[i_i_T];
					//for ( i_T = 0; i_T < SetT.Size; ++i_T )
					//{
						fMakeNewSeed_List.pElement[i_T] = true;
					}
					for ( i_sample = 0; i_sample < Amount_Sampling; ++i_sample )
					{
						for ( ; ; )
						{
							sampled_T = CMyRandom::DiscreteSample(vT_realP_T_Posterior.pElement, SetT.Size);
							if (fMakeNewSeed_List.pElement[sampled_T]) break;
						}
						Ideal_T_List.pElement[i_sample] = sampled_T;
					}
					fMakeNewSeed_List.Resize(Amount_Sampling);
					fMakeNewSeed_List.Fill(false);

					int RandomZ = Amount_Sampling / Ideal_Amount_Seed;
					int Amount_Seed = 0;
					for ( i_sample = 0; i_sample < Amount_Sampling; ++i_sample )
					{
						if (i_sample>0 && CMyRandom::NextInt(RandomZ)>0) continue;
						fMakeNewSeed_List.pElement[i_sample] = true;
						++Amount_Seed;
					}
					//sprintf_s(strTemp, "Seed Amount = %d;", Amount_Seed);
					//CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, true);
					//CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);
				}

				Ideal_T_Count = EPS;
				Selected_T_Count = EPS;
				for ( i_sample = 0; i_sample < Amount_Sampling; ++i_sample )
				{
					if (fMakeNewSeed_List.pElement[i_sample])
					{
						sampled_T = Ideal_T_List.pElement[i_sample];
					//if (i_sample==0) // || CMyRandom::NextInt((int)sqrt(1.0 * Amount_Sampling))==0) // Sample New T and New Base R
					//{
					//	sampled_T = CMyRandom::DiscreteSample(vT_realP_T_Posterior.pElement, SetT.Size);

						const CDistribution& tT_xyR_logP_Ra_given_T_RRa_FF_F = vT_xyR_logP_Ra_given_T_RRa_FF_F.pElement[sampled_T];
						for ( i_i_a = 0; i_i_a < (int)Waiting_Cell_List.Size; ++i_i_a )
						{
							i_a = Waiting_Cell_List.pElement[i_i_a];
						//for ( i_a = 0; i_a < IWH; ++i_a )
						//{
							int& t_Ra = oR_given_D_F.Result.pElement[i_a];
							const double* QRa_realP_Ra_given_T_RRa_FF_F = tT_xyR_logP_Ra_given_T_RRa_FF_F.pBase_realP(i_a);
							t_Ra = CMyRandom::DiscreteSample(QRa_realP_Ra_given_T_RRa_FF_F, IE);
						}

						// LY Modified 20090528
						//vT_logP_RR_R_T_given_FF_F.Clear(); vT_logP_RR_R_T_given_FF_F.PushBack(vT_logP_T_Prior);
						vT_logP_RR_R_T_given_FF_F = Extra_log_ppT_List;

						vT_CountZeroP_RR_R_T_given_FF_F.Fill(0);
						for ( i_i_a = 0; i_i_a < (int)Waiting_Cell_List.Size; ++i_i_a )
						{
							i_a = Waiting_Cell_List.pElement[i_i_a];
						//for ( i_a = 0; i_a < IWH; ++i_a )
						//{
							int t_Ra = oR_given_D_F.Result.pElement[i_a];

							for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
							{
								i_T = Waiting_T_List.pElement[i_i_T];
							//for ( i_T = 0; i_T < SetT.Size; ++i_T )
							//{
								double& logP_RR_R_T_given_FF_F = vT_logP_RR_R_T_given_FF_F.pElement[i_T];
								int& CountZeroP_RR_R_T_given_FF_F = vT_CountZeroP_RR_R_T_given_FF_F.pElement[i_T];
								const double* QRa_xyR_logP_RRa_Ra_given_T_FF_F = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T].pBase_logP(i_a);
								ttlp = QRa_xyR_logP_RRa_Ra_given_T_FF_F[t_Ra];
								if (ttlp > -MAXIMUM_DOUBLE_NUM_C)
									logP_RR_R_T_given_FF_F = CMyMath::log_Product(logP_RR_R_T_given_FF_F, ttlp);
								else
									++CountZeroP_RR_R_T_given_FF_F;
							} // i_T
						} // i_a

						sslp = -MAXIMUM_DOUBLE_NUM_V;
						for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
						{
							i_T = Waiting_T_List.pElement[i_i_T];
						//for ( i_T = 0; i_T < SetT.Size; ++i_T )
						//{
							if (vT_CountZeroP_RR_R_T_given_FF_F.pElement[i_T] == 0)
								sslp = CMyMath::log_Sum(sslp, vT_logP_RR_R_T_given_FF_F.pElement[i_T]);
						}

						CH_2 += sslp;
						to_divide += 1.0;
					}

					if (1) // Do Gibbs Sampling on R
					{
						int old_Ra;

						for ( i_i_a = 0; i_i_a < (int)Waiting_Cell_List.Size; ++i_i_a )
						{
							i_a = Waiting_Cell_List.pElement[i_i_a];
						//for ( i_a = 0; i_a < IWH; ++i_a )
						//{
							int& t_Ra = oR_given_D_F.Result.pElement[i_a];

#if 1
							old_Ra = t_Ra;
							{
								const CDistribution& tT_xyR_logP_Ra_given_T_RRa_FF_F = vT_xyR_logP_Ra_given_T_RRa_FF_F.pElement[sampled_T];
								const double* QRa_realP_Ra_given_T_RRa_FF_F = tT_xyR_logP_Ra_given_T_RRa_FF_F.pBase_realP(i_a);
								t_Ra = CMyRandom::DiscreteSample(QRa_realP_Ra_given_T_RRa_FF_F, IE);
							}

							int i_Effected_T;

							const int &Amount_Affected_T = AffectMatrix_Amount_T.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a];
							if (Amount_Affected_T < 0)
							{
								// Making AffectMatrix[i_a, old_Ra, t_Ra]
								int& ta_Amount_T = AffectMatrix_Amount_T.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a];
								CSimpleTypeArray<int>& ta_List_T = AffectMatrix_T.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a];
								CSimpleTypeArray<double>& ta_List_Delta_logP = AffectMatrix_Delta_logP.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a];
								CSimpleTypeArray<int>& ta_List_Delta_ZeroCount = AffectMatrix_Delta_ZeroCount.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a];
								ta_Amount_T = 0; ta_List_T.Clear(); ta_List_Delta_logP.Clear(); ta_List_Delta_ZeroCount.Clear();

								double t_logP;
								int t_ZeroCount;

								for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
								{
									i_T = Waiting_T_List.pElement[i_i_T];
									//for ( i_T = 0; i_T < SetT.Size; ++i_T )
									//{
									const double* QRa_xyR_logP_RRa_Ra_given_T_FF_F = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T].pBase_logP(i_a);
									t_logP = 0.0; t_ZeroCount = 0;

									ttlp = QRa_xyR_logP_RRa_Ra_given_T_FF_F[old_Ra];
									if (ttlp > -MAXIMUM_DOUBLE_NUM_C)
										t_logP = CMyMath::log_Product(t_logP, -ttlp);
									else
										--t_ZeroCount;

									ttlp = QRa_xyR_logP_RRa_Ra_given_T_FF_F[t_Ra];
									if (ttlp > -MAXIMUM_DOUBLE_NUM_C)
										t_logP = CMyMath::log_Product(t_logP, +ttlp);
									else
										++t_ZeroCount;

									if (t_ZeroCount==0 && abs(t_logP)<AffectThreshold_Ra)
										continue;

									++ta_Amount_T;
									ta_List_T.PushBack(i_T);
									ta_List_Delta_logP.PushBack(t_logP);
									ta_List_Delta_ZeroCount.PushBack(t_ZeroCount);
								}
							}
							const int* Qa_T = AffectMatrix_T.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a].pElement;
							const double* Qa_Delta_logP = AffectMatrix_Delta_logP.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a].pElement;
							const int* Qa_Delta_ZeroCount = AffectMatrix_Delta_ZeroCount.pElement[((size_t)(t_Ra + old_Ra * IE)) + AM_Size_IE_IE * (size_t)i_a].pElement;

							Ideal_T_Count += Waiting_T_List.Size;
							Selected_T_Count += Amount_Affected_T;
							for ( i_Effected_T = 0; i_Effected_T < Amount_Affected_T; ++i_Effected_T )
							{
								i_T = Qa_T[i_Effected_T];

								double& logP_RR_R_T_given_FF_F = vT_logP_RR_R_T_given_FF_F.pElement[i_T];
								int& CountZeroP_RR_R_T_given_FF_F = vT_CountZeroP_RR_R_T_given_FF_F.pElement[i_T];

								if (CountZeroP_RR_R_T_given_FF_F == 0)
									sslp = CMyMath::log_Minus(sslp, logP_RR_R_T_given_FF_F);
								logP_RR_R_T_given_FF_F = CMyMath::log_Product(logP_RR_R_T_given_FF_F, Qa_Delta_logP[i_Effected_T]);
								CountZeroP_RR_R_T_given_FF_F += Qa_Delta_ZeroCount[i_Effected_T];
								if (CountZeroP_RR_R_T_given_FF_F == 0)
									sslp = CMyMath::log_Sum(sslp, logP_RR_R_T_given_FF_F);
							}
#else
							old_Ra = t_Ra;
							//slp = -MAXIMUM_DOUBLE_NUM_V;
							for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
							{
								i_T = Waiting_T_List.pElement[i_i_T];
								//for ( i_T = 0; i_T < SetT.Size; ++i_T )
								//{
								double& logP_RR_R_T_given_FF_F = vT_logP_RR_R_T_given_FF_F.pElement[i_T];
								int& CountZeroP_RR_R_T_given_FF_F = vT_CountZeroP_RR_R_T_given_FF_F.pElement[i_T];

								const double* QRa_xyR_logP_RRa_Ra_given_T_FF_F = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T].pBase_logP(i_a);

								ttlp = QRa_xyR_logP_RRa_Ra_given_T_FF_F[old_Ra];
								if (ttlp > -MAXIMUM_DOUBLE_NUM_C)
									logP_RR_R_T_given_FF_F = CMyMath::log_Product(logP_RR_R_T_given_FF_F, -ttlp);
								else
									--CountZeroP_RR_R_T_given_FF_F;

								//tlp = CMyMath::log_Product(ttlp, vT_mat_logP_RRa_given_T_FF.pElement[i_T].pElement[i_a]);
								//vT_realP_T_given_RR_R_FF_F.pElement[i_T] = tlp;
								//slp = CMyMath::log_Sum(tlp, slp);
							}

							//for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
							//{
							//	i_T = Waiting_T_List.pElement[i_i_T];
							////for ( i_T = 0; i_T < SetT.Size; ++i_T )
							////{
							//	tlp = vT_realP_T_given_RR_R_FF_F.pElement[i_T];
							//	tlp = CMyMath::log_Product(tlp, -slp);
							//	vT_realP_T_given_RR_R_FF_F.pElement[i_T] = exp(tlp);
							//}
							//sampled_T = CMyRandom::DiscreteSample(vT_realP_T_given_RR_R_FF_F.pElement, SetT.Size);
							{
								const CDistribution& tT_xyR_logP_Ra_given_T_RRa_FF_F = vT_xyR_logP_Ra_given_T_RRa_FF_F.pElement[sampled_T];
								const double* QRa_realP_Ra_given_T_RRa_FF_F = tT_xyR_logP_Ra_given_T_RRa_FF_F.pBase_realP(i_a);
								t_Ra = CMyRandom::DiscreteSample(QRa_realP_Ra_given_T_RRa_FF_F, IE);
							}

							sslp = -MAXIMUM_DOUBLE_NUM_V;
							for ( i_i_T = 0; i_i_T < (int)Waiting_T_List.Size; ++i_i_T )
							{
								i_T = Waiting_T_List.pElement[i_i_T];
								//for ( i_T = 0; i_T < SetT.Size; ++i_T )
								//{
								double& logP_RR_R_T_given_FF_F = vT_logP_RR_R_T_given_FF_F.pElement[i_T];
								int& CountZeroP_RR_R_T_given_FF_F = vT_CountZeroP_RR_R_T_given_FF_F.pElement[i_T];

								const double* QRa_xyR_logP_RRa_Ra_given_T_FF_F = vT_xyR_logP_RRa_Ra_given_T_FF_F.pElement[i_T].pBase_logP(i_a);

								ttlp = QRa_xyR_logP_RRa_Ra_given_T_FF_F[t_Ra];
								if (ttlp > -MAXIMUM_DOUBLE_NUM_C)
									logP_RR_R_T_given_FF_F = CMyMath::log_Product(logP_RR_R_T_given_FF_F, +ttlp);
								else
									++CountZeroP_RR_R_T_given_FF_F;

								if (CountZeroP_RR_R_T_given_FF_F == 0)
									sslp = CMyMath::log_Sum(sslp, logP_RR_R_T_given_FF_F);
							} // i_T

#endif

							//vList.PushBack(CMyMath::my_log(-sslp));

							CH_2 += sslp;
							to_divide += 1.0;
						} // i_a

					} // Do Gibbs Sampling on R

					if (CMyRandom::NextInt(Amount_Sampling/100) == 0)
					{
						timCH_sfs.Stop();
						stc_2 = CExternalTools::GetTimeText(
							timCH_sfs.GetCumulativeTime() * (Amount_Sampling - i_sample - 1) / (i_sample + 1));
						timCH_sfs.Start();
						sprintf_s(strTemp, "F [%d/%d]; Progress [%02d%%], CH = %+8.4lf; TL = %s;",
							i_F, SetF.Size, 99*(i_sample+1)/Amount_Sampling, CH_1 + CH_2 / to_divide, stc_2.c_str());
						CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, false);
					}

					//if (1)
					//{
					//	int tt = i_sample + 1;
					//	while (tt>1 && !(tt&1)) tt >>= 1;
					//	if (tt==1 && i_sample+1>=2048)
					//	{
					//		sprintf_s(strTemp, "%6d\t%12.6lf", i_sample + 1, CH_2 / to_divide);
					//		CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);
					//	}
					//}
				} // i_sample

				timCH_sfs.Stop();

				CH_2 /= to_divide;

				//sprintf_s(strTemp, "%d  %.6lf", i_try, CH_2);
				//CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, true);
				////CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);

			} // Gibbs Sampling // or i_try
		} // Make CH_2
#endif

		timCH_2.Stop(); timCH_F.Stop();

		CH_F = CH_1 + CH_2;
		if (FL<0 || CH_F_L>CH_F)
		{
			CH_F_L = CH_F; FL = i_F;
			DstFixation = SetF.List.pElement[i_F];
		}
		if (FH<0 || CH_F_H<CH_F)
		{
			CH_F_H = CH_F; FH = i_F;
		}

		stc_1 = CExternalTools::GetTimeText(timCH_1.GetCumulativeTime()/(i_F+1));
		stc_2 = CExternalTools::GetTimeText(timCH_2.GetCumulativeTime()/(i_F+1));
		stc_F = CExternalTools::GetTimeText(timCH_F.GetCumulativeTime());
		sprintf_s(strTemp, "Done F [%d/%d]; Best F = <%+d,%+d> CH = %.4lf; TC = [%s];",
			i_F, SetF.Size, DstFixation.ob_x, DstFixation.ob_y, CH_F_L, stc_F.c_str());
		CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, false);

		sprintf_s(strTemp, "[%04d/%04d]\t<%03d,%03d>\t%-12.4lf\t%-12.4lf\t%-+8.4lf\t",
			i_F, SetF.Size, SetF.List.pElement[i_F].ob_x, SetF.List.pElement[i_F].ob_y, CH_1, CH_2, CH_F);
		CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, false);
		sprintf_s(strTemp, "%04d\t<%03d,%03d>\t%-+8.4lf\tTC = %s [ ( %s / %s ) per F ]; Per|a| = %.2lf%%; Per|T| = %.2lf%%;",
			FL, DstFixation.ob_x, DstFixation.ob_y, CH_F_L, stc_F.c_str(), stc_1.c_str(), stc_2.c_str(),
			100.0 * Waiting_Cell_List.Size / IWH, 100.0 * Selected_T_Count / Ideal_T_Count);
		CExternalTools::ConsoleTools::OutputTextToStream(oLog, strTemp, true);
	} // i_F

	stc_1 = CExternalTools::GetTimeText(timCH_1.GetCumulativeTime()/SetF.Size);
	stc_2 = CExternalTools::GetTimeText(timCH_2.GetCumulativeTime()/SetF.Size);
	stc_F = CExternalTools::GetTimeText(timCH_F.GetCumulativeTime());
	sprintf_s(strTemp, "Best F = %04d <%+04d,%+04d>, Expected H[T] ( CH[T|F] ) = %+8.4lf;",
		FL, DstFixation.ob_x, DstFixation.ob_y, CH_F_L);
	CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, true);
	oLog<<strTemp<<endl;
	sprintf_s(strTemp, "Time Cost = %s [ ( %s / %s ) per F ];",
		stc_F.c_str(), stc_1.c_str(), stc_2.c_str());
	CExternalTools::ConsoleTools::MakeCurrentLine(strTemp, true);
	oLog<<strTemp<<endl;

	// Output Distribution of CH[F]
	{
		CSimpleTypeArray<double> matCH;
		matCH.Resize(IWH);
		matCH.Fill(CH_F_H);
		int ix, iy;

#if 1
		CSimpleTypeArray<double> wl;
		wl.Resize(SetF.Size);
		for ( i_a = 0; i_a < IWH; ++i_a )
		{
			ix = i_a % IW;
			iy = i_a / IH;
			slp = -MAXIMUM_DOUBLE_NUM_V;
			for ( i_F = 0; i_F < SetF.Size; ++i_F )
			{
				tlp = wl.pElement[i_F] = 3.0 * CMyMath::my_log(func_e(
					ix - SetF.List.pElement[i_F].ob_x, iy - SetF.List.pElement[i_F].ob_y));
				slp = CMyMath::log_Sum(slp, tlp);
			}
			sslp = 0.0;
			for ( i_F = 0; i_F < SetF.Size; ++i_F )
			{
				tlp = wl.pElement[i_F] = CMyMath::log_Product(wl.pElement[i_F], -slp);
				sslp += exp(tlp) * DstList_CH_F.pElement[i_F];
			}
			matCH.pElement[i_a] = sslp;
		}
#else
		for ( i_F = 0; i_F < SetF.Size; ++i_F )
		{
			ix = SetF.List.pElement[i_F].ob_x;
			iy = SetF.List.pElement[i_F].ob_y;
			if (ix<0 || ix>=IW || iy<0 || iy>=IH) continue;
			matCH.pElement[ix + iy * IW] = DstList_CH_F.pElement[i_F];
		}
#endif
		i_sample = -1;
		for ( i_F = 0; i_F < SetF.Size; ++i_F )
		{
			ix = SetF.List.pElement[i_F].ob_x;
			iy = SetF.List.pElement[i_F].ob_y;
			if (ix<0 || ix>=IW || iy<0 || iy>=IH) continue;
			if (i_F == FL) i_sample = ix + iy * IW;
		}
		CSimpleTypeArray<int> mat_iCH;
		mat_iCH.Resize(IWH);
		for ( i_a = 0; i_a < IWH; ++i_a )
		{
			i_T = (int)(256.0 * (matCH.pElement[i_a]-CH_F_L) / (CH_F_H - CH_F_L));
			i_T = (i_T<0 ? 0 : (i_T>255 ? 255 : i_T));
			mat_iCH.pElement[i_a] = 255 - i_T;
		}
		int iiw = IW * Output_CW;
		int iih = IH * Output_CH;
		CImg<int> imgOutput;
		CImg<int> imgContext;
		imgContext.assign(iiw, iih, 1, 3);
		imgContext.fill(0);
		for ( iy = 0; iy < iih; ++iy ) for ( ix = 0; ix < iiw; ++ix )
		{
			i_F = ix/Output_CW + iy/Output_CH * IW;
			i_T = mat_iCH.pElement[i_F];
			imgContext.at(ix, iy, 0, 0) = i_T;
			imgContext.at(ix, iy, 0, 1) = i_T;
			imgContext.at(ix, iy, 0, 2) = i_T;
		}
		if (i_sample >= 0)
		{
			CRectBox<int> iBox;
			ix = i_sample % IW;
			iy = i_sample / IW;
			iBox.x = ix * Output_CW;
			iBox.y = iy * Output_CH;
			iBox.w = Output_CW;
			iBox.h = Output_CH;
			CExternalTools::DrawRectBox(imgContext, iBox, 255, 0, 0);
		}
		if (Answer_T_Index>=0 && Answer_T_Index<SetT.Size)
		{
			CRectBox<int> iBox;
			iBox.x = (ObjectTemplate.xLow + SetT.List.pElement[Answer_T_Index].base_x) * Output_CW;
			iBox.y = (ObjectTemplate.yLow + SetT.List.pElement[Answer_T_Index].base_y) * Output_CH;
			iBox.w = (ObjectTemplate.xHigh - ObjectTemplate.xLow + 1) * Output_CW;
			iBox.h = (ObjectTemplate.yHigh - ObjectTemplate.yLow + 1) * Output_CH;
			CExternalTools::DrawRectBox(imgContext, iBox, 255, 255, 0);
		}
		DrawGridImage(imgOutput, imgContext, IW, IH, Output_CW, Output_CH, Output_H_Extent);
		sprintf_s(strTemp, "Best F = %04d [%+04d,%+04d]; CH = %+8.4lf;",
			FL,	DstFixation.ob_x, DstFixation.ob_y, CH_F_L);
		CExternalTools::OutputTextToImage(imgOutput, strTemp, 4, imgOutput.height - Output_H_Extent + 2, 0, 0, 0, 255, 255, 255, 16);
		sprintf_s(strTemp, "%sf_given_R_F.jpg", strPN_Prediction_Output.c_str());
		CExternalTools::mySaveImage(strTemp, imgOutput);
	}

	DstBest_F = FL;
}

